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ABSTRACT 


The primary thrust of this thesis is the development of a computer-based 
simulation of an Orthogonal Frequency Division Multipexing (OFDM) based underwater 
acoustic communication system. The product will support the testing and evaluation of 
various digital signal processing algorithms applicable to underwater acoustic 
communication systems using OFDM as well as the study of the effects of the acoustic 
channel and communication system factors on the key parameters of the system such as 
bit error rate, received signal to noise ratio, frequency band of employment and overall 
system bit rate. The underwater acoustic sound channel is modeled using a physics based 
parabolic equation approximation. The simulation models the key components in the 
transmitter and receiver that contribute to the overall performance of the system. The 
results of the thesis provide expected values for system performance in terms of bit rate, 
bit error rate and received SNR for given frequency bands and are validated through 
comparison to theoretically derived expectations and to ocean testing of OFDM 
underwater communication systems. 
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EXECUTIVE SUMMARY 


Communication in the underwater battle spaee is essential to Netcentrie Warfare 
eoneepts in the 2U* eentury. Acoustie eommunieation holds promising eapabilities for 
eonneeting the undersea environment. The underwater aeoustie sound ehannel however, 
presents ehallenging obstaeles that must be overeome to allow satisfaetory aeoustie 
eommunieation. Possibly the most signifieant are the frequeney seleetive fading and 
multipath natures the underwater ehannel presents. 

Orthogonal Frequency Division Multiplexing (OFDM) is extremely robust to 
multipath frequeney seleetive fading ehannels and therefore is rapidly developing interest 
in application to underwater aeoustie communication. The primary thrust of this work is 
to develop a eomputer-based simulation of an OFDM based eommunieation system for 
the undersea environment. The result will then provide a foundation upon whieh the key 
parameters and eapabilities of underwater aeoustie eommunieation using OFDM ean be 
investigated. The immediate benefit of the work is the establishment of the aehievable 
eommunieation system parameters, eommon to all oommunieation systems, which are bit 
rate, bit error rate, bandwidth and the required signal to noise ratios. 

The simulation is performed in the Matlab environment where all major hardware 
in the transmitter and reeeiver are simulated. The underwater aeoustie sound channel is 
modeled using the physics based parabolie equation approximation. The model is the 
Monterey-Miami Parabolie Equation (MMPE) whieh employs the Split Step Eourier 
method. Use of a physios based model vioe the typioal stoohastio model for the aeoustie 
ehannel ensures a better realization of the phenomenon that is experienoed in the 
undersea environment. 

The simulation allows for any type of data to be transmitted from analog voioe 
signals to binary data. Signifieant flexibility is ooded into the simulation to allow for a 
wide range of parameter analysis and sensitivity testing. Reed Solomon ending and block 
interleaving are applied to aehieve several deeibels of forward error eorreetion ending 
gain. Quadrature amplitude modulation (QAM) is applied with eonstellation sizes of 8 to 
64 points. Double sideband suppressed earrier modulation is used to modulate the 
baseband signal to the desired transmission band. 

Bit rates achieved range from 4 to 10 kilobits per second (kbps) utilizing from 2 
to 6 kilohertz (kHz) of bandwidth depending on the size of the QAM eonstellation 
applied. The transmitted signals oeeupy frequeney bands between 4 and 14 kHz. 

The results are validated through a theoretical development of achievable 
performanee. In addition validation is made through eomparison to in water testing of 
OEDM. 

The results substantiate the feasibility of OEDM as an attraetive eandidate for 
underwater aeoustie communication. The product of the work provides validated system 
performanee expectations as well as a foundation upon whieh future work in this area ean 
be performed. 
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I. INTRODUCTION 


A. PURPOSE 

Netcentric Warfare has moved to the forefront of the Navy’s priorities for the 2U* 
century. Essential to the concept of Netcentric Warfare is the ability to communicate 
quickly, reliably, efficiently and in many cases covertly in all regions of the battle space. 
More specifically, key characteristics of the communication capability required are high 
data rates, low error rates, environmental robustness and minimally overt acoustic 
signatures. The undersea environment is arguably the most challenging region of the 
battle space in which to communicate effectively. Current and future goals for 
underwater wireless communication involve data flow between submarines, surface units, 
unmanned undersea vehicles (UUVs) and ocean surface and ocean bottom 
communication buoys in any combination via the underwater sound channel. Research 
and development in the fundamental communication methods and protocols that support 
Netcentric Warfare concepts will enhance the Navy’s ability to fight successfully in all 
areas of the battle space. 

Orthogonal Frequency Division Multiplexing (OFDM) possesses characteristics 
that have the potential to overcome the adverse effects of frequency selective multipath 
fading present in the underwater sound channel. Establishing groundwork in this area to 
address the feasibility of OFDM for underwater communication is essential to developing 
an underwater communication system that is beneficial to the Navy. 


B, GOALS 

The primary goal is to develop a computer simulation that will allow for research 
into the key parameters, capabilities, and limitations of OFDM in underwater acoustic 
communication. The simulation provides a ‘test bed’ for future study of other modulation 
methods, forward error correction coding, quantization, signal detection, decoding etc. In 


1 



developing and validating the simulation the feasibility of OFDM systems in underwater 
eommunieation is verified. There is one known experiment with OFDM in the 
underwater acoustic sound channel that is used as a standard against which the results of 
the simulation are validated. The work was done by Coatelan and Glavieux. [Ref. 1] 

While there has been considerable work using OFDM for communication all but 
the work of Coatelan and Glavieux [Ref. 1] and that of Kim and Lu [Ref. 2] has been 
outside the area of underwater acoustic communication. Kim and Lu [Ref. 2] present a 
new broadband underwater channel model that has the capability to model the Doppler 
time scaling effect as well as multi-path effects. The propagation path delay is modeled 
with two components that generate the simulated Doppler shift determined by the relative 
motion between the source and receiver. The simulation work by Kim and Lu [Ref 2] 
represents the only previous work in simulating an underwater acoustic communication 
system employing OFDM as the backbone. A physics based model is used to provide a 
realistic model of underwater acoustic sound channel. 

C. METHODOLOGY 

The simulation is run on a PC using Matlab code. The underwater sound channel 
is modeled using a physics based, parabolic equation (PE) model, namely the Monterey 
Miami Parabolic Equation model. 

The simulation is used to conduct experimentation to examine key parameters 
such as achievable bit rate and bit error rate, required SNR, transmission range, 
robustness to multipath propagation and frequency selective fading, and bandwidth 
requirements. The strengths and weaknesses of OEDM applied to underwater acoustic 
communications will be identified to assess its feasibility and merit for future work. 


D, BENEFITS 

Simulation is rapidly becoming a popular method for experimentation with 
technology-based systems. One of the biggest advantages of simulation is the ability to 
perform countless numbers of experiments under controlled conditions at a minimal cost. 
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Some of the biggest obstacles to simulation are the decision of what level of simulation is 
required to accurately capture the real world and validation of the simulation. With all 
simulations there must be an initial model. This model suffers from some simplifications 
that cause departure from reality. Later improvements will help the model to converge on 
the characteristics of the real world system that is being modeled. With that in mind many 
design decisions have been made in developing the simulation. The emphasis is on 
modeling the aspects of the system which have the greatest impact on the key parameters 
of interest. 

Finally, the benefit of this thesis is the ability to experiment with OFDM based 
communication systems such that the factors influencing the key parameters in 
communication can be understood. Understanding these factors will provide the 
knowledge necessary to build actual systems for underwater acoustic communication. 
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II. OFDM THEORY 


A, OFDM BASICS 

Modulating an orthogonal set of subcarriers with Q-ary complex QAM symbols 
generates an OFDM signal. The basic OFDM system is shown in Fig. 2.1 below. The 
input bit stream is buffered and converted to parallel complex QAM symbols in the QAM 
Modulator. 


Input Bit 
Stream^ R 

(bps) 


Output Bit 
Stream 



Figure 2.1. OFDM Basic Block Diagram. 


The subchannel, denoted as Xiy, where i is the subchannel index and k is the 


OFDM symbol index, can be independently modulated as in adaptive modulation 
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schemes or all subearriers may be modulated in the same manner. The parallel QAM 
symbols are then applied to the IFFT, whieh generates N parallel real samples. The real 
samples are then eonverted to serial format. In this simplified bloek diagram the reeeiver 
reverses the operations performed by the transmitter. 

The QAM modulator buffers the input bits into bloeks of size b, where b = RT,T 
is the symbol duration in seeonds and R is the overall bit rate in bits per seeond. The 
symbol rate then is Rs = 1/T symbols per second. The number of bits per subehannel 
determines the size of the QAM eonstellation, Q. For example with q = A bits per 
subehannel, Q = 2"^, whieh is 16 QAM. In OFDM the frequeney speetrum of the signal is 
ereated in the QAM Modulator and then the time domain signal results from passing the 
output through the IFFT operation. The eomplex QAM symbols derived from the binary 
eode words are sequentially plaeed in a subearrier band and OFDM symbol. Figure 2.2 
illustrates the organization of four OFDM blocks with eight subearriers. 
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Figure 2.2. Subcarrier and OFDM Symbol Arrangement. 


The IFFT provides an orthogonal transformation and maintains the energy 
contained in the input symbol [Ref 3] Written as an equation, this is 


1 N 2 

N 


n-l 


Due to the orthogonality of the IFFT the subcarriers each have an integer number 
of cycles within the FFT interval. The number of cycles in the FFT interval for a given 
subcarrier differs by one from that of adjacent subcarriers. Note that the terms subchannel 
and subcarrier are used interchangeably in OFDM terminology. Figure 2.3 illustrates the 
orthogonality of 3 subcarriers in an OFDM symbol. In this example each of the 
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subcarriers has the same magnitude and phase but this is not in general the case. As can 
be seen, each subcarrier has an integer number of eycles within the FFT interval. 



1. Subcarrier Orthogonality 

The real bandpass OFDM symbol ean be expressed as follows: 


2 / + 1 / 2 

x{t) = Re{ ^ - - fo))}, t,<t<t,+T 

‘ 2 

x(f) = 0,t <tQ t > t^+T 


( 2 . 2 ) 


where T is the symbol period, 2^ is the earrier frequency, N is the number of subearriers, 
di are the eomplex QAM symbols and to is the symbol starting time [Ref. 4]. From Eq. 
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2.2 it can be seen that the OFDM symbol consists of non-zero subcarriers over the 
symbol period T. As a result the frequency spectrum of the OFDM symbol consists of the 
convolution of dirac delta functions at frequencies corresponding to the individual 
subcarrier frequencies with the frequency spectrum of a unit amplitude square wave with 
duration T. The frequency spectrum of a unit amplitude square wave is smc{7rfr), which 
is zero at all frequencies that are a multiple of HT. As a result the spectrum of an OFDM 
symbol with overlapping sine functions is as shown in Figure 2.4. Note that at the 
maximum of each sine, i.e. frequency spectrum of each subcarrier, all other sine 
functions are zero. This is the point where the OFDM receiver calculates the spectral 
values for the individual subcarriers. This means that the OFDM receiver is able to 
demodulate each subcarrier free of interference from other carriers. Inter Carrier 
Interference (ICI) would exist and degrade the received signal if this were not the case. 
Note that this analysis is done in the frequency domain of the OFDM signal and therefore 
it is ICI rather than Inter Symbol Interference (ISI) that is avoided. [Ref. 4] 



Figure 2.4. Subcarrier Orthogonality via Subcarrier Spectra [after Ref 4]. 
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The real bandpass signal of Eq. 2.2 can also be expressed in complex baseband 


form as 



(2.3) 


The complex baseband signal is the inverse Fourier transform of N QAM 
symbols. Specifically each OFDM symbol is the inverse Fourier transform of the 
column s or blocks of Fig. 2.2. In practice the Inverse Fast Fourier Transform (IFFT) is 
applied to the QAM symbols to generate the discrete time OFDM signal written in 
equation form as 



(2.4) 


where N is the size of the IFFT and is equal to 2* A .[Ref. 2] 

Since the output of the IFFT must produce a real signal, the OFDM blocks of Fig. 
2.2 must be modified prior to application of the IFFT. From Fourier transform theory it is 
known that a complex frequency domain signal that is conjugate symmetric will produce 
a real time domain signal upon application of the inverse Fourier transform. Therefore 
since the N QAM symbols in each OFDM block represent the information to be 
transmitted, and are defined in the frequency domain, the blocks are conjugated and 
repeated to produce the conjugate symmetric counterpart of the original OFDM block. 
The result is an OFDM block of length N{=2* N ) complex QAM symbols. The first N 
QAM symbols represent the real frequency portion of the frequency spectrum and the 
remaining N QAM symbols represent the imaginary frequencies in the spectrum. Then 
upon application of the IFFT to each block of length N complex QAM symbols, a block 
of length N real time domain samples is generated which are converted to serial in 
preparation for transmission. 
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2. 


Guard Time and Cyclic Extension 


One of the major advantages of OFDM is its robustness to multipath propagation 
and the assoeiated delay spread that oeeurs. The delay spread is the time difference of 
arrival of transmitted signals that follow different paths from the source to the receiver 
and therefore have different propagation times. 

Dividing the input data stream into N subcarriers reduces the symbol duration by 

1/. To combat the effects of ISI a guard time is incorporated into each OFDM symbol. 
The duration of the guard time is chosen to be at least as long as the maximum expected 
delay spread. The maximum expected delay spread is the time difference of arrival 
between the first and last transmitted signals to arrive at the receiver. This ensures that 
multipath receptions of one symbol do not interfere with the receptions of the following 
symbol. 

One possibility for the guard time is to transmit no signal. However this would 
result in ICI, which is crosstalk between carriers that results from a loss of orthogonality. 
Recall that orthogonality means that all subcarriers have an integer number of cycles in 
the FFT period with adjacent subcarriers having a number of cycles that differ by one. 
Integration of a subcarrier which has a non-integer number of cycles in the FFT interval 
will result in a loss of orthogonality and therefore signal degradation due to ICI. This 
effect is illustrated in Figure 2.5 below. 

Figure 2.5 shows 2 subcarriers with the second subcarrier delayed. Note that the 
guard time consists of a zero signal. Also note that subcarrier one has exactly two cycles 
within the FFT interval while subcarrier two, due to its delay, does not contain an integer 
number of cycles within the FFT interval which will cause ICI as discussed above and 
indicated in Fig. 2.5. The shaded portion of subcarrier 2 indicates the region of the signal 
that causes ICI. One can see that if the subcarriers were extended cyclically in the guard 
interval that there would be an integer number of cycles in the FFT interval regardless of 
the multipath delay of subcarrier 2. The cyclic extension is implemented by simply 
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extending the length of each OFDM block following the IFFT operation by an amount 
greater than or equal to the maximum expected delay spread. 



Figure 2.6 illustrates the cyclic extension of the subcarriers in the guard interval. 
The benefit of using a cyclic extension in the guard interval is shown in Figure 2.7, which 
shows a 2 ray multipath environment with the amplitude of the second ray half the 
amplitude of the first ray. The second ray arrives after the first ray but its delay spread is 
less than the duration of the guard interval thus there is no ISI. 
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[after Ref 4]. 


In addition, it can be seen that there are an integer number of cycles for rays 
within the FFT integration time which maintains the orthogonality of the subcarriers and 
prevents ICI. The subcarriers in Fig. 2.7 are BPSK modulated. Therefore only two 
possible phases are possible and they differ by 180 degrees. The phase shifts that result 
are indicated in the figure. Note that the OFDM receiver would actually see the 
summation of all the received paths but the subcarriers are shown separately for 
illustration purposes. 
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environment [after Ref. 4], 


As discussed above, the maximum delay spread must not exceed the guard 
interval incorporated in the OFDM signal or ISI and a loss of orthogonality will result 
causing ICI. Van Nee and Prasad [Ref. 4] show that with the delay spread greater than the 
guard interval by less than about three percent, the signal interference is acceptable but at 
around 10% delay spread excess the interference is unacceptable for communication 
purposes. However OFDM is still significantly more robust to delay spreads that exceed 
the value for which the system was designed relative to that for single carrier systems as 
error propagation causes abrupt performance degradation in single carrier systems [Ref. 

4]. 
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B. 


CHANNEL ANALYSIS 


In general communication systems, the input is considered to be the transmitted 
signal and the output is the received signal. The received signal is the transmitted signal 
modified by the channel and corrupted by additive noise. The channel is usually band 
limited and assumed to be linear. It may be either time varying or time invariant. This 
concept is illustrated in Figure 2.8 below, where x(n) is the input signal, y(n) is the 
output signal, h(n) is the channel impulse response and w(n) is additive gaussian noise. 



In OFDM, each of the N subchannels is transmitted with its own carrier. Consider 
an OFDM symbol defined as 



(2.5) 


where k refers to the OFDM block or symbol. As discussed above, Xk must be modified 
to be conjugate symmetric prior to applying the IFFT. Therefore Xk becomes 
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( 2 . 6 ) 


where A^= 2* N . The IFFT then modulates the N subearriers resulting in the real time 
domain signal Xk, 



(2.7) 


The channel is assumed to have a finite impulse response of length Z. Therefore a 
cyclic extension of length L is applied to the OFDM symbol xt. The OFDM signal 
becomes 






V 

cyclic extension 




IFFT{X,) 





Transmitted OFDM Signal 


( 2 . 8 ) 


The output from the channel is the convolution of the channel impulse response 
with the input signal. Written as an equation this is 
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(2.9) 



[n] 


[n)® h{n) 


where 0 represents cireular convolution and yk(n) is the received OFDM block. It is 
well known that for continuous time signals the linear convolution in the time domain is 
equivalent to multiplication in the frequency domain. In order for this to hold in discrete 
time one of two conditions must be met. Either the block length N must be infinite or at 
least one of the convolved sequences must be periodic with period N. Observe that the 
input OFDM signal is periodic with period N due to the application of the cyclic 
extension. Therefore the received frequency domain signal is the product of the input 
OFDM frequency domain signal and the frequency response of the channel. In equation 
form the kth received OFDM symbol, with the subscript k omitted for simplicity, is 


Y^HX 


i i 


for i = 1 , 2 . 



( 2 . 10 ) 


where the subscript i refers to the subcarrier. Note that the use of the cyclic prefix reduces 
the achievable data rate by a factor of ^ . However for most applications N»L so 

the loss is negligible relative to the identified benefits of using the cyclic prefix. Most 
OFDM systems use the cyclic prefix [Ref. 5 and 6]. 


The received signal is 





Received OFDM Signal 


( 2 . 11 ) 
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As indicated in Eq. 2.11 the last N values oiyk are the eonvolution of the OFDM 
input signal with the ehannel impulse response. Therefore the received QAM symbols in 
the k*’’ bloek are written as Eq. 2.12 where the last N QAM symbols are the eomplex 
eonjugate of the first N QAM symbols and are disearded in the reeeiver. 



and h(n) is zero padded to length N in taking the IFFT. As a result the effeet on the 
ehannel is that it is also divided into narrow subehannels eorresponding to the subearriers. 

For large N , the subehannels are very narrow and there is negligible change in H^ and 
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versus frequency. In other words the frequency selective fading that may be 

present is eliminated. Subchannel independence can be shown provided the additive noise 
is Gaussian(whether the noise is colored or not) [Ref 3]. This result is illustrated in Fig. 
2.9 below. 

The rectangles represent the subchannels that are used by the subcarrier bands. 
Each subchannel is assumed flat for large N. In Fig. 2.9, only 11 subcarriers exist 
however in practice the number of subcarriers is typically between 500 and 2000. With 
this in mind one can see that the frequency response within each of the subchannels 
would be flat as the width of the subchannel would be drastically reduced. 
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III. OFDM SIMULATION MODEL 


The development of a eomputer-based simulation of an underwater OFDM 
eommunieation system is the primary thrust of this thesis. The model simulates most of 
the significant hardware that would exist in an actual system. This chapter describes the 
details of the model and the theory behind the methods applied. The intent is not to 
exhaustively model every component in an actual system but to model those components 
and features of a real system that contribute the dominating characteristics of the 
system’s performance. It is also understood to be unfeasible and foolish to attempt to 
develop a model from scratch that will have the ability to experiment with every possible 
design and operating feature of the system. With that in mind, the model was developed 
to be comprehensive enough to explore some of the basic specifications common to any 
communication system immediately but then to also provide a test bed upon which 
further modeling and experimentation could take place. 


A, METHODOLOGY 

The development of each block in the model followed essentially the same path. 
First the topic of interest was explored in the literature to the extent required. The extent 
of literature research was a function of the number of choices available to accomplish the 
task, the experience of this author in the subject and the relative importance of the block 
to the entire model. For example, the modulation and demodulation blocks were only 
mildly researched while the topic of forward error correction was researched much more 
extensively. Figure 3.1 provides a block diagram of the entire simulation. 
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B, 


TRANSMITTER 


The transmitter covers all aspects of signal processing and transmission from 
information signal development to Quadrature Amplitude Modulation (QAM), to Double 
Side Band Suppressed Carrier (DSBSC) modulation of the signal. Each of the transmitter 
blocks of Fig. 3.1 is discussed in the following numbered subsections. 


1, OFDM Parameter Determination 


In designing a communication system there are several parameters of interest that 
must be determined. The final design will be a compromise between conflicting 
requirements to achieve the optimal performance. Parameter determination is performed 
in the simulation by the Matlab function parameters.m contained in Appendix C. In 
general there are three main parameters to determine: bit rate, R, channel delay spread, 
tds, and bandwidth, W. Delay spread is a function of the channel and therefore represents 
the opening argument in designing the system. The OFDM system must be designed such 
that the tolerable delay spread of the system is greater than the delay spread of the 
channel to avoid signal degradation. Therefore the guard time, tg, is typically two to four 
times the maximum expected delay spread. That is. 



(3.1) 


The delay spread for the underwater acoustic sound channel is a function of but 
not limited to bottom type, transmission range, transmitter and receiver depths and 
frequency band. Coatelan and Glavieux [Ref 1] report measured delay spreads in Table 
3.1 below. The table identifies the transmission range between the receiver and the 
transmitter, the seafloor composition and the associated tds measured. 
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Bottom 

Transmission 

Delay Spread 

Type 

Range (m) 

(msee) 

Rocky 

< 1000 

< 10 

Sludgy 

1500-2500 

<10 

Sandy 

1500-2500 

>20 

Unknown 

1430 

~20 


Table 3.1. Typical Delay Spread Values. 


The next step is to determine the OFDM symbol duration. The OFDM symbol 
duration is the sum of the cyclic prefix duration and the FFT interval duration. Since the 
signal in the cyclic prefix does not contain any new information it is desirable to have the 
FFT interval duration much greater than the guard time to maximize the achievable bit 

rate. However there is a limit. As T increases, the subcarrier spacing, , decreases for 

a given bandwidth. Smaller subcarrier spacing results in greater sensitivity to phase errors 
and frequency offset. In practice the symbol duration is typically 5 or more times tg, 
therefore 


T > ; 

5 * t 


g 


For channels with frequency selective fading, smaller values are preferred 
as this ensures the subchannel frequency response is fiat. Conversely, for channels with 
large Doppler spreading characteristics larger A/^^ values are better. 

Figure 3.2 below presents correlator output results measured by Coatelan and 
Glavieux [Ref. 1]. The results in Fig. 3.2 correspond to the 1430 m acoustic channel of 
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Table 3.1 with a sludgy seafloor and a depth of 90 m. Note that the delay spread in Fig. 
3.2 is approximately 20 msee and the first path reeeived is not neeessarily the most 
energetie one [Ref. 1]. 



With the guard time and symbol duration established, the number of subearriers 
N can be determined. There are two approaches that can be taken. The first is starting 
with a known bandwidth constraint. In this case the number of subcarriers is simply the 
available bandwidth, W, divided by the subcarrier spacing, therefore 


¥sc \ 

The other approach, which is the one used in the simulation, is to start with a 
desired bit rate, R. Then the number of subcarriers is simply R divided by the number of 
bits transmitted by each subcarrier, bsc- That is. 
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(3.4) 



a. Example Calculation 


Suppose that a system is to be designed that would provide a 10 kbps data 
rate and the maximum expected delay spread for the acoustic channel of interest is 27 
msec. That is, 

R = 10 kbps 
= 27 msec 

let = 4 ( 27 ) = 108 msec 

let T = 6t^ = 6 ( 1 08) = 648 msec 


Assume that the system is capable of transmitting 6 bits per subcarrier, 
which corresponds to 64-QAM, then 


N = 



\0E?> 

6 


= 1667 


subcarriers 


However since the QAM symbols on each of the subcarriers will be 
modulated using the IFFT, it is computationally more efficient to have the number of 

subcarriers equal to a power of two. Therefore N is set to 2048 subcarriers but only 1667 
of the subcarriers have non-zero values with the remaining subcarriers providing 
oversampling to avoid aliasing. In essence the signal itself is unaffected and changing N 
only affects the internal signal processing efficiency. Next the FFT interval duration is 
calculated as. 


RfFT ~ ^ —108 — 540 msec 
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Another requirement is that there be an integer number of samples within 
both the FFT interval, Tfft, and the entire OFDM symbol, T, in order to maintain 
orthogonality of the subearriers. For this example there must be 2048 samples within the 
FFT interval. The FFT sampling rate, Fs_fft, is then. 


F 


N 2048 


s,FFT 


FFT 


540 


= 3.793 kHz 


However the number of samples within the OFDM symbol interval, T, is, 

* r = 3.793 * 540 = 2457.6 


whieh is not an integer. Therefore the subearriers will not be orthogonal and signal 
degradation will oecur. The remedy is to adjust the sampling rate slightly. The sampling 
rate is increased to achieve an integer number of samples within the FFT interval and the 
OFDM symbol interval. For this example becomes 3842.6 Hz, which yields 

Up = F^ * J" = 3842.6 * 648F — 3 = 2490 samples 
npFT = F^ ppp * Tppp = 3842.6 * 540F — 3 = 2075 samples 
—f^s,FFT^tg = 3842.6* 108F— 3 = 415 samples 


where rij ,, Uppj, and are the number of samples in the OFDM symbol interval, the FFT 

interval and the guard interval respectively. Note that the OFDM symbol is oversampled 
by increasing the sampling rate since there are 2075 vice 2048 samples within the FFT 
interval. Synchronization, which will be discussed in the receiver section, will determine 
which of the 2075 samples in the FFT interval are to be used for signal reconstruction. As 
a result of changing the sampling rate. 


T = = 0.5330 sec 

^ s,FFT 


=-^ = 0.115 sec 

^s,FFT 
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and finally, 


¥sc =t;^t = 1-876 Hz 

W = N*Af^^= 2.0264 kHz 
2, Signal Generation 

In underwater communication several possible data forms may be transmitted 
including video, voice and text format messages. For the simulation, voice data in .wav 
format is used as the information signal. Voice data was chosen as it requires an analog 
signal to be converted to binary format via sampling, quantization and decimal to binary 
conversion. Once the data is in binary form the original form of the data is irrelevant to 
the OFDM system. Figure 3.3 presents the information signal used in the simulation. The 
upper plot is the discrete time signal. The sampling rate of the analog voice signal is 
11.025 kHz and the data is quantized using 8 bits per sample. The lower plot is the 
frequency spectrum of the transmitted signal. 
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Figure 3.3. Information Signal. 
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3. 


Quantization 


The sampled information signal is then quantized using uniform quantization. The 
level of quantization is adjustable but 8-bit quantization is primarily used in the 
simulation whieh results in 256 quantization levels. Quantization is performed by the 
Matlab funetion quantization.m contained in Appendix C. 

The output of the quantizer is a stream of decimals from 0 to M-1 where M is 2*”'® 
and bits is the number of bits per binary digits per quantization level. The decimal value 
corresponds to the quantization level for which the sampled value of the information 
signal was nearest. The quantizer also outputs a codebook that identifies the quantization 
level corresponding to the decimal value for signal reconstruction. Figure 3.4 clarifies the 
quantization process with a simple 2-bit example. 
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Figure 3.4. Quantization Illustration. 
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4. 


Decimal to Binary Conversion 


Decimal to binary conversion is performed by the Matlab code 
binary conversion.m contained in Appendix C. The process involves straightforward 
conversion of a base 10 integer to the corresponding base 2 representation. 

5. QAM Buffer 

QAM modulation involves the mapping of bits to a QAM constellation. The size 
of the constellation, Q, is a function of the number of bits per QAM symbol. 

g = 2^ P.5, 

In general the number of bits per QAM symbol will not be the same as the 
number of bits per quantization level. Therefore it is necessary to buffer the bits prior to 
mapping to a QAM constellation. The Matlab code qam buffer.m contained in Appendix 
C performs the buffering operation. This involves converting the bit stream from the 
decimal to binary conversion process into a matrix of bits. Each row in the matrix 
represents a QAM symbol in binary form and the number of columns per row is q. The 
number of rows is a function of the size of the data to be transmitted. Since the matrix 
must be square in order to convert all bits to a QAM symbol, up to q-\ zeroes are added 
to the data. 


6, Reed Solomon Encoder 

Due to the multipath, frequency selective fading nature of the underwater acoustic 
sound channel the subcarriers will arrive with different amplitudes. Some subcarriers may 
experience severe enough fading that they are lost entirely. As a result the bit error for 
those subchannels will approach 14 and will dominate the BER performance of the entire 
system[Ref. 4]. 

Eorward error correction coding is used to combat the severe fading affect. The 
coding is applied across all subcarriers such that the errors of the weak subcarriers can be 
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corrected. The performance of the system is then a function of the average received 
power vice the received power of the weak subcarriers [Ref 2], 

In the simulation block coding is applied, specifically Reed Solomon block 
coding which is performed at the symbol level. Coding involves the addition of redundant 
symbols to increase the Hamming distance between symbols. In general a message of 
length k is added thereby creating a code word of length n. The coding rate r is defined as 
k/n. The minimum Hamming distance, dmin, is the minimum number of different symbols 
between any 2 code words. 

For a given code with an associated dmin, t errors can be corrected by the code 
where t is 



where the floor function denotes rounding down to the nearest integer[Ref 2]. The 
minimum Hamming distance has an upper bound defined by the number of redundant 
symbols as 

(3.7) 

In terms of bit level codes, only repetition and single parity check codes reach this 
upper bound. Reed Solomon is a class of non-binary codes that reach this upper 
bound. [Ref 4] 

Reed Solomon codes are defined for code words of length n symbols associated 
with m bits per symbol where 

(3.8) 

The required message length, k, is related to message and the number of bits per 
symbol as 

(3.9) 
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Associated with the number of symbol errors the eode ean eorreet is the number 
of bits the eode is eapable of correeting, tt, defined as 


< m * floor 


^ n - 


(3.10) 


Equation 3.10 is true provided the bit errors oecur within the maximum 
eorreetable number of symbol errors. For example eonsider a eode designed to be able to 
eorreet 3 symbol errors with 4 bits per symbol. The eode will not be able to eorreet an 
arbitrary oeeurrenee of 4 bit errors as they may oecur in 4 different symbols. As a result 
Reed Solomon eodes are useful for ‘bursty’ error ehannels such as the multipath, 
frequeney selective fading underwater aeoustie channel sinee the errors are eoneentrated 
in a few of the subearriers that are deeply faded. 


The Reed Solomon eneoder reeeives QAM eoded deeimals in serial format from 
the QAM buffer. The vector of decimal symbols must be reordered to a A: eolumn square 
matrix. Therefore subsymbols may need to be added to ereate a square matrix, in whieh 
case they are added to the beginning of the veetor. 


7. Block Interleaver 

The frequeney seleetive fading eauses varying amplitudes of the subearriers 
where some of the earriers will be unreliable due to severe fading. As a result burst errors 
viee random errors will oeeur over the weak subearriers. Most EEC eodes are designed 
for random errors not burst errors[Ref. 4]. The transmitter permutes the data sueh that the 
adjaeent symbols are separated by several bloeks after interleaving. The reeeiver 
performs the inverse operation. 

As diseussed Reed Solomon eodes are designed to eorreet a eertain number of 
errors per bloek. The interleaving must distribute the subsymbols aeross several bloeks 
sueh that the bursts eaused by the severely faded subearriers appear random in the 
reeeiver after de-interleaving. Therefore rather than having bursts of symbol errors within 


32 




each block corresponding to the severely faded subcarriers, the errors within each block 
will be randomly distributed which is optimal for the Reed Solomon decoding algorithm. 

The interleaving process is essentially a reordering of the data and is performed 
by the reshape.m command in Matlab. 

The input data to the interleaver is an n column matrix where the data is ordered 
in time sequence across the rows. Therefore the sequential order of the data in the matrix 
begins at row 1, column 1, then proceeds across row 1 to column n, then to row 2, 
column 1 etc. 

The reshape function takes the data column wise such that the subsymbols are 
interleaved in frequency. The process is illustrated in Fig. 3.5 with n equal to 7 and k 
equal to 3. In the simulation the length of the input vector A is on the order lO"^, n is at 
most 64, and the number of active subcarriers per OFDM symbol is typically less than 
10 . Therefore the interleaving is over several blocks as desired. 
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Figure 3.5. Reed Solomon Eneoder and Bloek Interleaver. 


8, QAM Modulator 

QAM Modulation maps the QAM coded decimals representing points on a QAM 
constellation into real and imaginary components. The real component is called the in- 
phase component, and the imaginary component is called the Quadrature component of 
the signal. The mapping from decimal to real and imaginary components is illustrated in 
Fig. 3.6 below. 
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Figure 3.6. 16 QAM Constellation with In-phase and Quadrature Mapping. 


The Cartesian coordinates of in-phase and Quadrature components are converted 
to polar coordinates prior to creating the OFDM symbol blocks. The QAM Modulator 
creates the frequency spectrum of the signal using the complex QAM symbols now in 
magnitude and phase components. The goal is to generate blocks of complex QAM 
symbols such that applying them to the IFFT operation results in a real signal. The IFFT 
is given by 



(3.11) 


where A is the length of the DFT. The discrete time signal x„ will be real provided the 
following condition of conjugate symmetry holds for the frequency domain signal A, 
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(3.12) 



The OFDM blocks are then formed by sequentially filling the active subcarriers in 
each block with the complex QAM symbols. 

In the parameter determination algorithm the number of active subcarriers, Nact, 

and the number of zero subcarriers, Nzer , are determined. The zero subcarriers are 
required to prevent aliasing from signal power near the Nyquist sampling rate which 

corresponds to N . The number of data blocks, Ntiock required to transmit a given amount 
of data is given by the number of complex QAM symbols to be transmitted divided by 
the number of active subcarriers rounded up to the nearest integer. As a result of 
rounding up to get the Nuock, there is a surplus of active subcarriers relative to the number 
of QAM symbols to be transmitted. In future work this surplus could be filled with 
protocol information. The surplus subcarriers are filled with zero symbols in the 
simulation. 


The original design placed the required number of zeroes at the beginning of the 
first block. However, this was found to be a poor choice as it results in an uneven power 
distribution across the frequency spectrum. The current design spreads the zero symbols 
across the frequency band used for transmission. 


to 


The number and location of zero symbols to be inserted are determined according 



(3.13) 


where is the total number of zero symbols to be inserted and ^ is the total 

number of complex QAM symbols to be transmitted. In general will not be a multiple 

of Nbiock, therefore the number of zero symbols inserted will not be the same for each 
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OFDM block. The following method is used to determine the number of zeroes inserted 
into eaeh OFDM bloek. 


f 


^ col = floor 


N 




0 


V ^ block J 


'^col,\ ^0 block 1) * ^co/ 


(3.14) 


where j is the number of zero symbols inserted into the first OFDM block and Zcoi is 

the number of zero symbols inserted into the remainder of the OFDM blocks. The 
process is illustrated in Fig. 3.7 below. 

Applieation of Eq. 3.12 to the N QAM symbols in eaeh of the OFDM blocks 
results in OFDM symbols whieh are eonjugate symmetric. The conjugate symmetry will 
ensure that each of the OFDM bloeks of complex QAM symbols will result in N{ = 2N ) 
real samples upon applieation of the IFFT. The result is shown in Fig. 3.8 with the same 

parameters as introduced in Fig. 3.7. The symbol ^ in Fig. 3.8 refers to eomplex 
eonjugation. 
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Figure 3.7. QAM Modulator. 


The QAM Modulator also creates the pilot blocks which are used for channel 
estimation and equalization. Although the theory behind channel estimation and 
equalization will be discussed later, the creation of the pilot blocks is presented here. 

In simple terms the pilot blocks are symbols known to the receiver that are used to 
estimate the channel frequency response. In the simulation two pilot blocks are 
implemented and make up the first two OFDM symbols that are transmitted. The pilot 
blocks are made up of pseudo random complex OFDM symbols which are determined as 
follows. 

The data used to create the pilot blocks comes from the data blocks illustrated in 
Fig. 3.7. Each complex QAM symbol is drawn from the data blocks and placed in the 
first pilot block as follows, 
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Xp,i{c,r)= X,Xr,c) for 


r — N act “1“ 1 5 N act 5 ... 1 

< 

c = N N 1 •( 

L'^ block ? block 

3.15) 

Equation 3.12 is applied to the first N eomplex QAM symbols ereating a length 
Abloek of complex QAM symbols that is eonjugate symmetric. Reversing the order of 
the first pilot bloek ereates the seeond pilot bloek. 

The final effeet of the QAM Modulator is that by adding the eomplex eonjugate 
eounterpart to each of the N eomplex QAM symbol bloeks the input data is effeetively 
interpolated. Typieally interpolation involves the insertion of one or more zeroes in a 
time sequence with the effeet of inereasing the sampling rate and ‘eompressing’ the 
frequeney speetrum of the sequenee. In the QAM Modulator the interpolation eonstant is 
2 as the number of data points is doubled. The result is that the sampling rate, Fs,fft, is 
doubled. Therefore the new sampling rate in the EFT interval is, 

(3.16) 
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9. 


IFFT/OFDM Modulator 


The IFFT is used to modulate the complex QAM symbols to the subcarrier 
frequencies. In other words each subcarrier is modulated by the corresponding complex 
QAM symbol. 



The IFFT can also be expressed in matrix form as 



(3.17) 


(3.18) 


(3.19) 


where e and e (jnxm ^]^gj.gfQj.g ^ j^nxm conjugate symmetric 

as defined in Eq. 3.12. 


Note that C refers to the complex subspace and R refers to the real subspace and 
the M refers to the number of columns in the matrix. In the simulation M is Ntiock- The 
output from the IFFT is a matrix of real samples with each column representing a 
baseband real OFDM symbol. Columns 1 and 2 are the pilot symbols and the remaining 
column s contain the information to be transmitted. 


41 





10, Cyclic Extension 


The need for the eyclie extension was diseussed in Chapter II. The 
implementation of the eyclie prefix in the simulation is quite simple. First the required 
length of the cyclic extension is determined by 

L = round(t^ 

where L is the number of discrete time samples that must be prefixed to each real OFDM 
symbol. 

Recall that from the IFFT is an x Nbhck matrix of real time samples 

therefore the cyclic extension results in a matrix that is N+L x Nbhck with the first L 
samples being identical to the last L samples in each column. Therefore the resulting 
matrix is 



The columns of the matrix of Eq. 3.21 represent the real baseband OFDM signals 
that are to be transmitted. 
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11. Parallel to Serial Converter 


The parallel to serial eonverter simply reshapes to a veetor of length 

{Nbiock+'2)*{N+L). In the simulation, round-off error causes the output from the IFFT to 
contain small imaginary components. Therefore the real part of the serial time samples is 
taken. 


12. Double Side Band Modulator 


After serial conversion the signal is baseband and must be modulated to the 
desired transmission band. The modulation is performed in discrete time to take 
advantage of the ability to use high order digital filters. 


In order to prevent aliasing the signal must be interpolated. The interpolation 
factor, /, is determined as 



(3.22) 


where fc is the carrier frequency. 


Interpolation is performed by inserting 7-1 zeroes between each sample. The 
interpolated signal can be expressed as 



where is the baseband signal andyfm) is the baseband signal following interpolation 
by the factor 7. The increased sampling rate is then 


F 


5, mod 



(3.24) 
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Figure 3.9 illustrates the effect of interpolation on the signal. The signal of Fig. 3.9 is for 
a 5 kbps transmission rate using 16-QAM. There are 1024 subcarriers with 750 active 
subcarriers which results in a bandwidth, W, of 1.501 kHz. The baseband signal is 
interpolated by a factor of 4 as shown in the lower graph of Fig. 3.9. 


Baseband OFDM Signal 



Frequency (Hz) 



Figure 3.9. Interpolated Baseband OFDM Signal. 


The interpolated signal, y(m), must be filtered to remove the unwanted spectral 
energy above the OFDM signal bandwidth, W, as seen in Fig. 3.9. A digital Finite 
Impulse Response (FIR) filter is used. The FIR filter is designed using a Hamming 
window with the following parameters. 
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I <»„ =-= pass band digital irequency 
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mod 


1.375(2;rlf) , 

1 0), = -= stop band digital frequency 




F 

5, mod 

8;r 

0). -co„ 


= filter order (number of coefficients is +1) 


The filter parameters are shown for a general FIR filter designed using the 
Hamming window method in Fig. 3.10. The variable CO, is the digital frequency defined 
as 



(3.25) 


where/refers to the analog frequency and F^ is the sampling rate used to sample the 
continuous time signal. The output of the FIR filter is defined as, 

yf{n)=h{n)^y{n) ^ 

where h(n) is the impulse response of the FIR filter and symbol * refers to linear 
convolution. 


Figure 3.11 shows the effect of the FIR filter on the interpolated baseband 
spectrum shown in the lower graph of Fig. 3.9. The interpolated baseband frequency 
spectrum is repeated in the upper graph of Fig. 3.11. The lower graph shows the 
frequency spectrum of the interpolated baseband frequency spectrum following filtering. 
The transmission parameters are the same as for Fig. 3.9. The frequency response of the 
FIR filter is shown in Fig. 3.12. 

After interpolation and filtering, the signal is ready for modulation to the desired 
carrier frequency. The modulation is performed by discrete time domain multiplication 
with a cosine modulating signal. The modulating signal is defined as 
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(3.27) 



COS 


^ 27tf^ 

F 

\ s, mod J 


where n is defined from 1 to the number of samples contained in the signal to be 
modulated, yf(n). 



Figure 3.10. General FIR Filter (Hamming Window Design). 
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The frequency spectrum of the cosine modulating signal consists of dirac delta 
functions at the absolute value of the carrier frequency on both the real and the imaginary 
frequency axis. The modulation is expressed in discrete time as 


Tmod («)="*(«)>< T/(«) 


(3.28) 


Baseband Interpolated OFDM Signal 
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Figure 3.11. Interpolated Baseband OFDM Signal following FIR Filter. 


The effect of modulation ony/'is shown in Fig. 3.13. The upper graph in Fig. 3.13 
repeats the frequency spectrum of the baseband OFDM signal prior to interpolation. The 
center graph is the frequency spectrum of the modulating signal which illustrates the 
dirac delta functions at the carrier frequency. Finally the lower graph is the frequency 
spectrum of the modulated OFDM signal. This is the form of the signal that is transmitted 
through the channel. In practice the signal would be applied to a Digital to Analog 
converter (DAC) and an associated low pass filter (LPF). For the simulation these 

47 




























components are assumed to be a part of the channel transfer function. Similarly, an 
Analog to Digital Converter (ADC) and associated LPF would be used in practice to 
return the signal to discrete time following transmission through the channel. These are 
also assumed to be included in the channel transfer function. 
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Figure 3.12. FIR Filter Frequency Response. 
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Figure 3.13. Modulation Signals. 


The DSB Modulator completes the transmitter for the simulation. The next effect 
on the signal in practice would be the channel. However in the simulation the channel 
transfer function is not applied until after demodulation in the receiver for reasons that 
will be discussed later. 

C. MMPE CHANNEL MODEL 

MMPE is the Monterey-Miami Parabolic Equation model. It provides an 
approximate solution to the Helmholtz equation based on the efficient split-step Eourier 
(SSE) marching algorithm. The Helmholtz equation is derived from the solution of the 
linear, three-dimensional, lossless, homogeneous wave equation describing the 
propagation of small amplitude acoustic signals in an ideal fluid [Ref. 6]. The velocity 
potential is assumed to have a time-harmonic solution, or in other words, the acoustic 
field is generated by a continuous wave source. The value of a time harmonic field at any 
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point and time in space is a function of only one frequency. Time harmonic solutions are 
important because solutions to the wave equation with arbitrary time dependence can be 
found through application of Fourier transform techniques to the time harmonic solutions 
[Ref. 7]. 


The following development follows the work by Houdeshell. [Ref. 8] The linear 
wave equation is given in Eq. 3.29 with the complex acoustic pressure field given as Eq. 
3.30, 


V + 


1 d^p 


0 



(3.29) 


(3.30) 


where P is the complex acoustic pressure field and p is the amplitude of the acoustic 
pressure. The Helmholtz equation is then 



(3.31) 


A cylindrical coordinate system is used in the solution. Therefore the position 

vector, f , is a function of range, r, depth, z, and the azimuthal angle, ^ . The azimuthal 
dependence of the acoustic field is minor and ignored in the solution. Therefore p is a 
function of range and depth only, i.e. p(r,z). 


Eet the complex acoustic pressure be defined as 


p 1 

[r ,z ,(Ot) 

1 = p \ 

(r , z )e ^ ^ ^ 


(3.32) 


and 



(3.33) 
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where the / /— term accounts for azimuthal spreading and u(r,z) defines the two 
/ Vr 

dimensional acoustic pressure field. Substituting Eqs. 3.32 and 3.33 into the Helmholtz 
equation yields the uncoupled azimuthal approximation, which is 


(3.34) 



The following operators, denoted by the subscript, op, are defined as 

(3.35) 


P = 

°p / d r 


and 


Qop - + ^ + i 


(3.36) 


where 


£ — n 


{331) 


and 



(3.38) 


Then Eq. 3.32 can be rewritten as 



If u is defined as 



(3.40) 


(3.39) 


then the outgoing wave can be shown to satisfy [Ref 9] 
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where the propagator, $(r), which is a unitary operator, steps the solution out in range. 
The split-step Fourier method is used to apply the unitary operator. The operator. Hop, is 
separated by using the wide-angle approximation [Ref. 10]of the square root operation 
such that 
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Therefore the propagator may be expressed as [Ref 11] 



The PE/SSF method applies the operator 



in the kz domain where the operator, , is defined as 


II 

1 

\ 



2 
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^ ko j 



The following convention is used for the FFT 


(3.47) 


(3.48) 


(3.49) 


(3.50) 


(3.51) 


(3.52) 
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(3.53) 


y>(^)= FFT {^(k.)} 

and 

yf(k, )= IFFT V(z)} 


Finally, the PE/SSF method is implemented as 


At / 

/ \ ~J\—U \z,r+Ar) 

y/{z,r+Ar)=e ^ FFT< 


-A^Ujz,r) , J 

e 2 yAz,r)Vf 



JJ 


(3.54) 


The MMPE model provides multiple outputs to the user. The primary output for 
use with the OEDM simulation is the eomplex acoustic pressure. The MMPE model used 
in the simulation is a 2D model in range and depth since the azimuthal dependence was 
neglected thereby omitting the third dimension. The solution is a grid of points from 
which the complex acoustic pressure can be evaluated. The complex acoustic pressure 
output for a given point in the grid is a function of frequency. The application of this data 
to the OEDM signal is discussed in the receiver section covering the application of the 
acoustic channel. 

The source depth for the shallow water (100 m) acoustic channels was 30 m and 
100 m for the deep water (340 m) acoustic channels. The receiver depth varied with the 
results of the acoustic channel evaluation. The receiver depth was chosen by determining 
the depth with the greatest average acoustic pressure amplitude across the frequency band 
evaluated. The sound speed profde (SSP) used for all acoustic channels is from the 
Atlantic Ocean just off the eastern coast and is shown in Eigure 3.14. 
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Figure 3.14. Sound Speed Profde. 


D. RECEIVER 

The funetional blocks in the receiver were presented in Fig. 3.1. For the most part 
the blocks of the receiver simply undo the work of the transmitter to convert the 
information back to its original form. Additional blocks in the receiver include the 
channel estimation and equalization blocks, additive white gaussian noise (AWGN), 
synchronization and bit error analysis. 

1. AWGN Addition 

The simulation assumes the background noise present in the channel is AWGN. 
The signal power is calculated as 
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(3.55) 



i = 1 _ 

K 


where K is the length of the reeeived signal. The noise power to be added is then 
ealeulated as, 



(3.56) 


where SNR is the desired SNR in dB. The reeeived signal is then 



( ) 

1 = x\ 

( ) 

l + w( 

[n) 


(3.57) 


where x(n) is the noise free signal at the input to the receiver, which is the output of the 
transmitter and w(n) is AWGN. 


2, DSB Demodulator 

The received signal is first demodulated using a cosine signal, as presented in Eq. 
3.26, then the demodulated signal is filtered and decimated. The EPF is designed such 
that spectral components above the known carrier frequency plus the bandwidth of the 
transmitted signal are removed. Then decimation by D is performed to reduce the number 
of samples and the sampling rate. 

The reader should note that in the current flowchart of the simulation, some 
processing of the received signal occurs before the acoustic channel model is applied. 
This is because without source or receiver motion the MMPE model can be applied 
directly to the complex QAM symbols in each of the OFDM blocks. In reality the 
channel effects would alter the signal after transmission and before reception. In future 

work source and receiver motion will be modeled. At that point the application of the 
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MMPE model will take plaee between transmission and reception. Therefore the blocks 
such as the DSB Modulator and Demodulator have been left in the simulation even 
though they are not essential in the current form. 

Signal detection theory is not addressed in this thesis. As such it is assumed that 
the signal is detected by some means and the stream of data containing the signal is 
processed. The simulation is designed such that the received signal of finite length 
contains the transmitted signal of some length less than the length of the received signal. 
The information is assumed to be received in an AWGN environment. The DSB 
Demodulator demodulates the entire received signal to baseband. The synchronizer is 
responsible for identifying where the information lies within the received signal. Only the 
portion of the received signal containing information, i.e. the transmitted signal, is 
processed in the remainder of the receiver. 


The demodulating signal is expressed as 



(3.58) 


where n is defined from 1 to the length of the received signal, Xr(n). The demodulated 
signal is then the product of the received signal with the demodulating signal, that is. 



( \ 

yn] 

1 = xj 

{ \ 
yn] 

\y.d{n) 


(3.59) 


The signal is baseband and real at this point. Next the signal must be decimated to 
reduce the number of samples and the sampling rate by the decimation factor, D. First the 
signal must be passed through an antialiasing LPF to remove spectral components above 

. The FPF is again designed using the Hamming window method with the following 
parameters: 
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I CD„ =-= pass band digital irequency 




s. mod 


1.25(2;rlf) , 

|<y„ =-= stop band digital irequency 




p 

i,mod 

S;^ 

C0-(0„ 


= filter order (number of coefficients is Nj^ +1) 


where referenee ean be made to Fig. 3.10 for the general filter eharaeteristies. The output 
of the filter is 



[n] 

1 = h{n] 


[n) 


(3.60) 


where h(n) is the impulse response of the LPF and * represents linear eonvolution. The 
final step in the DSB Demodulator is deeimation of the signal. Deeimation by a faetor D 
results in the expression 



(3.61) 


The effeets of the DSB Demodulator are illustrated in Fig. 3.15. The signal 
parameters of Fig. 3.15 are the same as those of Fig. 3.9. The bit rate is 5 kbps using 16 
QAM and the deeimation faetor is D is 4. The baseband bandwidth is 1.501 kHz. 
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Figure 3.15. DSB Demodulator Signal Proeessing Effeets. 


3, Synchronizer 

In order to demodulate the subearriers an OFDM receiver must first perform 
synchronization. There are two tasks related to synchronization. The first is determining 
the location of the OFDM symbol blocks. The second is estimating and correcting for 
frequency offset of the subcarriers. In addition coherent receivers also require the carrier 
phase of the DSB Modulator to be synchronized. In the simulation the phase of the DSB 
Modulator is assumed to be synchronized by the OFDM receiver. 

Recall from the discussion on OFDM theory that the subcarriers will remain 
orthogonal in the receiver only if the transmitter and receiver use exactly the same 
frequencies. If not then ICI will occur. Practical oscillators possess phase noise, which 
results in varying subcarrier center frequencies since the frequency is the time derivative 
of the phase. In practical systems with non-ideal oscillators this is an unavoidable source 

of ICI in OFDM systems. However Van Nee and Prasad [Ref. 4] indicate that this source 
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can be minimized to aceeptable levels through a well-designed system. In the simulation 
phase noise is not modeled and therefore assumed negligible. The sensitivity of OFDM 
systems to phase noise and frequeney offset is one of the major disadvantages of OFDM 
relative to single carrier systems. In single earrier systems phase noise and frequency 
offset degrade SNR but do not produce interference [Ref 4]. Van Nee and Prasad [Ref 
4] provide a detailed discussion of these two effects in OFDM systems. 

Frequency offset is eaused not only by phase noise of local oscillators in the 
transmitter but also by Doppler shift due to non-zero relative motion between the source 
and the receiver. The simulation assumes the transmitter and receiver are stationary. 
Therefore frequeney offsets are not addressed in this thesis. 

While OFDM systems are relatively sensitive to frequency issues they are quite 
robust to timing errors. Symbol timing may vary by as much as the duration of the guard 
interval, tg, without producing ISI or ICI. However an optimal timing instant exists which 
provides for maximum multipath robustness. Any deviation from the optimal timing 
instant reduces the delay spread tolerance the system was designed to handle. OFDM 
systems attempt to minimize the timing error relative to the guard interval. 

The cyelic extension is used to perform symbol synchronization in the simulation. 
Recall that the first tg samples reeeived will be the same as the last tg samples for eaeh 
OFDM block due to the cyelie extension. Therefore correlation streams of data tg samples 
long with streams of the same length but delayed by the duration of the FFT interval will 
identify the loeation of the OFDM bloeks. 

The synehronization bloek diagram is presented in [Ref. 4] and modified slightly 
for the purposes of the simulation to the form shown in Fig. 3.16. 


60 



c(«) 



m = l,2...2L-l 


Figure 3.16. Synchronization Flowchart [after Ref 4]. 


Recall that tg is the guard time in continuous time and L is the discrete time 
number of samples corresponding to tg for the applicable sampling rate. Figure 3.17 
illustrates the structure of an OFDM signal in serial form to clarify the cross correlation 
used in the synchronizer. Fig. 3.17 identifies the observation intervals of length N+L 
which slide over the received signal in discrete time. The first and last L samples in the 
observation interval will be the same when the interval falls exactly across an OFDM 
symbol labeled as frames in Fig. 3.17. At this point the cross correlation terms will be the 
greatest. When the observation interval includes samples from two different frames such 
that the first and last L samples are independent the eross correlation funetion will have a 
peak much lower than when the interval is exactly over one frame. In fact in the limit as 
the number of samples over which the eross eorrelation is performed is very large the 
ratio of the sidelobe-to-peak amplitude will go to zero [Ref 4]. Figure 3.18 provides a 
realization of the correlation signal identified in Fig. 3.16 for an OFDM signal with six 
OFDM blocks where the first two blocks are the channel estimation blocks. 
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Figure 3.17. OFDM Signal Structure with Synchronizer Observation Intervals. 



The peak finding algorithm uses five threshold levels to determine the existence 
of peaks in the correlation signal. The threshold levels are adjustable and based on the 
maximum value found in the correlation signal. In the realization of Fig. 3.18 the 
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threshold levels are 0.4, 0.45, 0.50, 0.55 and 0.60 of the maximum of the correlation 
signal. The algorithm requires that for a value to qualify as a peak it must first define an 
inflection point in the correlation signal. Second the value must register as a peak at all 
threshold levels. The impulse response and frequency response of the acoustic channel 
model applied to the transmitted signal in realizing the correlation signal in Fig. 3.18 is 
shown in Fig. 3.19. The impulse response of Fig. 3.19 clearly indicates the existence of 
three distinct receptions due to the multipath environment. Also note that the frequency 
response of Fig. 3.19 illustrates frequency selective fading and nearly linear phase 
response. 


Acoustic Channel 


Frequency Response 





Figure 3.19. Acoustic Channel Filter Plots for Realization of Fig. 3.17. 


Van Nee and Prasad [Ref 4] discuss optimal timing by analyzing the correlation 
signal power. The idea is to take the samples that lie within the OFDM symbol that also 
have the greatest power and therefore the greatest SNR. In the simulation the signal 
power is not analyzed since the correlator synchronizer, while computationally intensive, 
was very successful at identifying the boundaries of the OFDM symbols. In most cases 
the correct symbol boundaries were found to within a few percent error relative to the 
length of the cyclic prefix, L. Future work might look at the benefits of using the optimal 
timing instant. Additionally, the correlator synchronizer used here is computationally 

intensive and a less demanding algorithm would be desirable for practical application. In 
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that case the optimal timing instant would be of greater value as a less intense timing 
synchronization method would likely be less aceurate. A low eomplexity 
synchronization method that seems to be of value in OFDM is presented by van de Beek 
et al.[Ref. 12], The synehronization is performed by the Matlab funetions peakfinder.m 
and synehronizer.m eontained in Appendix C. 

4, Composite Description of Serial to Parallel Conversion through 
FFT/OFDM Demodulation 

The serial to parallel (S/P) eonverter takes the redueed length signal from the 
synehronizer and eonverts it to parallel form. The resulting matrix of samples will have 
the number of rows eorresponding to the bloek length determined by the synehronizer 
and the number of eolumns will be the number of OFDM symbols found by the 
synehronizer. Reeall that eaeh peak in the eorrelation signal eorresponds to the start of an 
OFDM symbol. The matrix of real time samples from the S/P converter eontain the eyclie 
prefix whieh is removed from eaeh OFDM symbol, i.e. eolurnn in the matrix, by the 
eyclie extension removal bloek. 

The OFDM Demodulator simply performs an N point FFT of the matrix from the 
eyelie extension removal bloek. The FFT ean be expressed as 

(3.62) 


(3.63) 

where g (jnxn ^ j^nxm ^ ^nxm resulting matrix of 

eomplex samples represent the transmitted complex QAM symbols. The OFDM symbols 
making up eaeh eolurnn of the matrix are eonjugate symmetrie and only the aetive 
subearriers are of interest for reeeiving the information in the signal. 



The FFT ean also be expressed in matrix form as 
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5. 


Acoustic Channel Application 


The output of the MMPE model used in the simulation is in the form of a eomplex 
transfer function. The transfer function is applied directly to the complex QAM symbols 
present after the OFDM demodulator. The generic output from the MMPE model is a 
grid of points in depth and range. The desired output is chosen for a given point. The 
complex transfer function represents the complex acoustic pressure as a function of 
frequency present at the chosen point in the grid. 

The transfer function contains values for the frequencies between fc and fc + 

■ This band represents the frequencies that would lie in the real frequency region 

of the baseband signal corresponding to the signal modulated to fc. Recall that only the 
active subcarriers carry non-zero signal power and therefore determine the bandwidth of 
the transmitted signal. In order to apply the transfer function to the complex QAM 
symbols it must be made conjugate symmetric by application of Eq. 3.12. The complex 
QAM symbols after application of the transfer function is 


Y ^ X .* H 


(3.64) 


where Xg and tfe therefore Yg The variable ^ 

represents the matrix of complex QAM symbols after the OFDM demodulator, 

represents the complex transfer function and Y- is the output. The .* in Eq. 3.64 

Y 

indicates array multiplication such that each column of ^ is weighted by the vector 

H . The number of column s in XandY is + 2 due to the 2 blocks of complex 

QAM symbols used for channel estimation and equalization. 


6. Channel Equalization 

The channel equalization scheme simply attempts to invert the effects of the 

channel such that the complex QAM symbols after equalization will be the same as those 
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that were generated in the transmitter. The first two OFDM symbols are used to estimate 
the channel transfer function. The generation of the channel estimation blocks was 
discussed in the QAM Modulator section. The channel estimate is 


(3.65) 

where A- represents the received complex QAM symbols in the channel estimation 
blocks, X is the complex QAM symbols present in the QAM Modulator of the transmitter 
and Hest is an A by 2 matrix of complex values representing estimates of the magnitude 
and phase of the channel transfer function at frequencies corresponding to the subcarrier 
center frequencies. Array division is implied in Eq. 3.65. The values ofXare assumed to 
be known at the receiver. Each row of Hest is averaged to give a vector of values 
representing the estimate of the channel transfer function. That is, 

(3.66) 

In the simulation the estimates of the channel transfer function are performed over 
the active subcarriers only rather than over all N as using the zero subcarriers would 
result in division by zero in Hest- 

Einally the complex QAM symbols are equalized by defining 

(3.67) 

where the division is array division such that the columns of X- are weighted by the 
inverse of . 

A realization of the channel estimates is shown in Eig. 3.20. The magnitude of the 
channel estimate for the two channel estimation blocks and the mean of the two estimates 
is plotted versus the EET bin number which is analogous with subcarrier number or 
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frequency. The plot is discontinuous at high and low bin numbers due to division by zero 
that would result where the zero symbols are interspersed to evenly shape the power 
distribution as discussed in the QAM Modulator section. The plots are over the active 
subcarriers only. Note that Hmean in Fig. 3.20 refers to the mean of the two channel 
estimates Hsys(l) and Hsys(2). 



Figure 3.20. Channel Estimates, A=2048, 16 QAM, 5 kbps,/c = 6 kHz, Nact = 750. 


7, QAM Demodulator 

The QAM demodulator transforms the received complex QAM symbols to points 
on a g-ary QAM constellation. Therefore the output will be integers from 0 to Q-\. 

The first step is to strip the first N act + 1 symbols from each block as these are 
the symbols that carry the transmitted information. This implies that the receiver knows 
the number of active subcarriers. This information is assumed to be contained in protocol 
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information that would be carried as overhead in the transmission system. Recall that the 
zero symbols that were added in the QAM modulator are for data management purposes 
and could just as well be protocol information. Protocol is not addressed in this thesis but 
is assumed to be present in making some assumptions that are required in the receiver. 

The zero symbols are removed and the magnitude and phase of each symbol is 
determined and converted to Cartesian coordinates. The real and imaginary components 
of the received symbols correspond to the in-phase and Quadrature components 
respectively. The operation of the QAM Demodulator up to this point is illustrated in Fig. 
3.21. 



Figure 3.21. QAM Demodulator. 


The next step is to map the in-phase and Quadrature components to Q-ary QAM 

constellation points. The mapping of constellation points to in-phase and Quadrature 

components was explained in the QAM modulator. In that case the constellation point 
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was exact and resulted in integer values for the in-phase and Quadrature components. 
However in mapping back to QAM constellation points the in-phase and Quadrature 
components will not in general be integer values. Therefore an algorithm is required to 
determine the constellation point to which the received symbol most closely corresponds. 

This operation is performed by the Matlab function qaskdecomod.m contained in 
Appendix C. A 2 by g matrix containing the in-phase and Quadrature values 
corresponding to each point in a ^-ary QAM constellation is created. The in-phase and 
Quadrature components of the received symbols are compared to the vectors in the 
constellation matrix to determine the constellation point that is nearest to the received 
symbol. The received symbol is mapped to this constellation point by its corresponding 
integer value. The process is illustrated in Fig. 3.22. The example in Fig. 3.22 is for a 
constellation size of 4 for simplicity although 4-ary QAM constellations are not used for 
data transmission in the simulation. The error analysis is performed for each received 
symbol resulting in a mapping of an integer from 0 to Q-\ for each in-phase and 
Quadrature component pair. Therefore the output of the QAM Demodulator is a vector of 

integer values of length A act * ^block ^0 ■ 
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8, Composite Description of Block De-Interleaving through Error 
Analysis 

The Block De-interleaver simply undoes the interleaving of data performed by the 
Block Interleaver in the transmitter. It is assumed that the receiver knows the code word 
length of the Reed Solomon code as the code vector supplied to the deinterleaver is 
reshaped into a matrix with n columns in preparation for Reed Solomon decoding. 

The Reed Solomon decoder applies an algorithm in the Matlab communications 
toolbox, rsdeco.m, in decoding. Recall the code word length is n and the message length 
is k. Therefore the input matrix has n column s and the output matrix has k column s where 
each row corresponds to a codeword when coded and a message when decoded. 

The integers representing QAM symbols must be converted to binary to continue 
the signal reconstruction. Each integer is converted to a binary word with q bits where q 
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is the number of bits per QAM symbol. The output of this block is a vector of binary 
digits. 

The quantization buffer simply reshapes the bit stream to a matrix with eight 
columns, since 8-bit quantization was performed in the transmitter. As stated earlier the 
level of quantization is adjustable and determined by the variable bits in the transmitter. 
In general the bit stream is reshaped to a bits column matrix in preparation for signal 
reconstruction. 

Signal reconstruction consists of two steps. The first is the conversion of the 
binary words of length bits to integers. Second the integer values are mapped to the 
corresponding value of the information signal in the transmitter. Refer back to Fig. 3.4 
which is the illustration of the quantization process. 

The BER is determined by comparing the transmitted and received bit streams. 
An error occurs when the received bit is different than the corresponding transmitted bit. 
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IV. RESULTS 


The results of the thesis are divided into 2 major seetions. The first section 
develops the theoretical performance of an OFDM system. The second presents the 
results of the OFDM simulation for several combinations of system and acoustic channels 
parameters. 

A, SYSTEM THEORETICAL PERFORMANCE 

The theoretical performance of the system is first analyzed by considering a single 
carrier system. An OFDM system is equivalent to a set of independent and ISI free single 
carrier QAM systems. Therefore a large part of the analysis of an OFDM system can be 
performed by analyzing a single carrier QAM system. The single carrier analysis is 
extended to the multicarrier, OFDM, system to complete the analysis. 


1. Single Carrier Analysis 


by 


The energy per two-dimensional symbol in a square QAM constellation is given 



(4.1) 


where d is the distance between points in the QAM constellation and Q is the size of the 
QAM constellation. All points in the constellation are equally likely and the constellation 
is centered at the origin and therefore has zero mean. Equation 4.1 is exact only for 
square QAM constellations but is still a good approximation for the average transmit 
energy for non-square QAM constellations [Ref. 1]. For an ISI free channel with gain, 
\h\ , the probability of two dimensional symbol error is approximately 



(4.2) 
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where dmin is the minimum distance between QAM constellation points at the output of 
the channel and Qer/x) is the error function, each of which are defined as [Ref. 3] 



and 


(4.3) 



(4.4) 


The probability of symbol error per dimension, which is PJ2, for communication 
systems should be less than 10'"^. Therefore the following requirement must be met. 



n.%(dB)+rSdB)-rAdB) 


(4.5) 


The variable is the margin, which provides for additional performance in order 
to overcome unpredictable channel degradations and y^ is the coding gain. Note that 
when y„ = yc=^ the system is uncoded and has no margin. Therefore 11.8 dB is the 

argument required in the Q(.) function to get = 10”"^. Coding increases dmin and 

therefore reduces the 11.8 dB by the amount of coding gain. Using a margin increases the 
11.8 dB required by the amount of margin desired. [Ref. 3] 

Equation 4.1 can be rewritten as 


1 

h 

o^ 

H 

2 

mm 


(4.6) 


which leads to the definition of SNR gap, or normalized SNR, Y , 
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(4.7) 

where is the variance or power of the AWGN. 

Therefore, for the target of = 10“"^, the following expression can be written 
for the SNR gap, 

?>Y{dB)=u.%+r„-r^{dB)\ ,, 3 , 

Then by taking the base 2 log of Eq. 4.6 and substituting into Equation 4.7 for 
dmin, the number of bits that can be carried by QAM at the target Pel2 is calculated as 

bp^ = logjie) = logji 1 + — I ,,,, 

where the channel output SNR is defined as 

(4.10) 

Note that Eq. 4.9 is the Shannon-Hartley theoretical channel capacity with the 
SNR reduced by a factor of the SNR gap. As the SNR gap goes to 1 (0 dB) the achievable 
data rate of the system approaches capacity. As can be seen, the SNR gap is a measure of 
the degradation of the system from optimal. The number of bits that can be carried by a 
QAM system as computed using Eq. 4.9 will not be an integer but is typically rounded 
down to the nearest integer. If a square QAM constellation is desired then it is rounded 
down to the nearest even integer. The SNR gap for any QAM system with probability of 
one dimensional symbol error equal to 10'"^ is 

r = 7.029 + 7m “ 7c (4.11) 
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2. 


OFDM (Multi-Carrier) Analysis 


With the single earrier analysis complete, multi-carrier analysis can now be 
performed as an extension of the single carrier results. In a generic multicarrier system 
the probability of error is an average of the probability of error for each of the 
subchannels. Therefore the overall probability of error will be dominated by the weakest 
subcarriers. [Ref 3] 


This weakness was pointed out earlier in the thesis. Forward error correction 
coding with block interleaving is used to average out the probability of error over the 
active subcarriers. Therefore in the following analysis the probability of error is assumed 
equal over all subchannels. Directly applying the single carrier analysis yields the 
following results for the i* subchannel 


1 

1 

3 ' 

1 

1 

H > 

2 ^ 2 
d, 

^ j- — ^ — 

4c7 ? 

Aaf 


(4.12) 


where the subscript i denotes any values that are subchannel dependent. The maximum 
number of bits that can be carried on the subchannel with a given margin and coding 
gain is 



(4.13) 


where 



(4.14) 


The subsymbol energy, , is held constant over all subchannels that are active 

and is zero on the inactive subchannels. The water pouring distribution (see Proakis [Ref. 

13] or Gallager [Ref. 14]), is a better energy distribution but Cioffi has determined that 
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the on/off energy distribution is very close to optimal yet is much easier to compute and 
implement. [Ref. 3] 

Figure 4.1 plots curves of bpe versus received channel SNR for three different 
values of one-dimensional probability of symbol error, Pei. Figure 4.1 assumes that the 
coding gain is equal to the margin which may be zero. Therefore actual results in 
simulation or that may be seen in real world application may shift left or right as a 
function of the margin and coding gain. If the coding gain is greater than the margin the 
curves will shift left thereby predicting greater numbers of bits per subcarrier at a given 
received SNR. The coding gain achieved in the simulation is not calculated but is 
assumed to be around 6 to 8 dB which will be shown to be accurate in the following 
results. 



Figure 4.1. Channel Bit per Subcarrier Capacity. 


The bit rate, R, of the total system is then simply the product of the active 
subcarriers with the number of bits per subcarrier divided by the OFDM symbol period. 
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T. Since all active subcarriers have the same number of bits, bsc, which is rounded down 
to the nearest integer below bp j , the system bit rate is 



(4.15) 


Figure 4.2 plots system bandwidth versus bit rate for several sizes of QAM 
constellations. 



Figure 4.2. Bandwidth vs Bit Rate for Q-ary QAM OFDM Signals. 


B. SIMULATION RESULTS 

The remainder of this chapter presents the results of the simulation over the 
acoustic channel models developed. The main parameters of the acoustic channel are the 
water depth, the bottom roughness in meters of standard deviation, the range between the 
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receiver and the transmitter and the frequency band over which the model is evaluated. 
The water depths evaluated are either 100 m or 340 m. Three different bottom roughness 
values of 0, 2 and 4 m standard deviation were evaluated. The ranges evaluated were 
either 2 km or 4 km. 

Two of the available forms of output from the MMPE model are included for 
every acoustic model evaluated. The first is a plot of transmission loss as a function of 
frequency and depth. The second is a plot of transmission loss versus depth and time. 
These plots show the ray paths present in the water column over time such that for a 
given depth the multipath arrivals are apparent. 

While each acoustic model was evaluated for ranges of 2 and 4 km, only the plots 
for 2 km ranges are provided. Figure 4.3 shows the transmission loss for the acoustic 
channel that is 100 m deep and has a bottom roughness of 4 m. The effect of the bottom 
roughness on the acoustic energy is clearly seen in Figure 4.3 and can be compared with 
the transmission loss curves of the same acoustic channel with less bottom roughness, 
which are Figures B.l and B.3, to further illustrate the effect. Note the interference 
patterns and blending of transmission near the interface that is much less pronounced as 
the bottom roughness goes down. Figure 4.4 shows the acoustic energy arrival as a 
function of time and depth for the same acoustic channel. The severe multipath 
environment present in shallow water acoustic channels is readily apparent in Figure 4.4. 
All other acoustic channel plots are found in Appendix B. 

The simulation was performed for several combinations of bit rate, center 
frequency and constellation sizes. Each combination was simulated 10 times and then the 
mean of the result was taken. The primary output of the simulation and this thesis are 
plots of BER versus received SNR. Recall that the simulation is performed with the 

margin, , set equal to zero. Accepted values for 7. range from 6 to 12 dB, with 
Cioffi [Ref 3] recommending 6 dB. 
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Figure 4.3. Transmission Loss, Bottom Roughness 4 m, Water Depth 100 m. 


Transmission Loss (dB re Im) 



The margin can be considered as an adjustment for simulations to shift the results 

towards what can be expected in reality. Recall also that the coding gain, 7c > is estimated 

to be 4 to 6 dB. Therefore, since the coding gain is included in the results of the 
simulation, in order to check the results versus the theoretical prediction of Fig. 4.1 the 
BER versus received SNR must be shifted to the right by approximately 11 dB. 
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Figure 4.4. Acoustic Energy Arrival, Bottom Roughness 4 m, Water Depth 100 m. 


The raw results are provided in Appendix A. In these plots each data point 
indicates the result of one run of the simulation. The solid red line is the mean of the 10 
runs performed. The composite results are provided in Figures 4.5 through 4.8 for QAM 
constellation sizes of 8, 16, 32 and 64. Each solid line in the figures is the mean of the 
corresponding series of 10 runs plotted individually in Appendix A. Figure 4.5 presents 
the results for all evaluations using 8-ary QAM. 
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Figure 4.5. 8 QAM Composite Results. 



Figure 4.6. 16 QAM Composite Results. 
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Figure 4.7. 32 QAM Composite Results. 



Figure 4.8. 64 QAM Composite Results. 
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Figure 4.9 provides an appreciation for the nature of the decoded signal and its 
frequency spectrum with respect to the information signal and its frequency spectrum. 


Transmitted Information Signal 



Figure 4.9. 16 QAM Signal Realization,= 6 kFlz, R = 5 kbps, SNR = 5.3 dB 

BER= 1.6 E-2. 


The BER is unacceptable by communication system standards. The information 
signal for the simulation is a voice stream. Playback of the received signal at the BER of 
Eigure 4.9 contains minor static but is clearly understandable. The zoomed-in subplots 
show the deviation of the received information signal spectrum from the transmitted 
information signal spectrum. One can see that where there is significant power in the 
frequency spectrum of the transmitted information signal the received spectrum is very 
near that of the transmitted spectrum. Comparing the time domain plots of the transmitted 
and received information signals one can see the deviations of the received signal relative 
to the original information signal. Eigure 4.10 shows the received QAM constellation for 


84 





















































































the realization of Figure 4.9 after channel estimation and equalization but prior to Reed 
Solomon decoding. 


Received QAM symbols prior to Reed Solomon decoding 
16 QAM, BER=1.6E-2 
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Figure 4.10. 16 QAM Signal Realization,yj; = 6 kFiz, R = 5 kbps, SNR = 5.3 dB 


BER= 1.6 E-2. 


Figures 4.11 and 4.12 provide a contrast to the realizations of Figures 4.9 and 4.10 
at a greater received SNR of 7.4 dB and the resulting smaller BER of 2.5 E-4. The BER 
is still greater than the acceptable limit of 1 E-4. However one can see that the received 
signal is essentially an exact replica of the transmitted information signal. Also the 
powerful effect of the EEC can be seen from the noisy constellations of Figures 4.10 and 
4.12. 
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BER = 2.5 E-4. 


Note that the zoomed plots of the frequency spectrums of Eig. 4.11 show that the 
received spectra is essentially the same as that of the transmitted information signal 
spectrum even in the areas of the spectrum where the transmitted information signal has 
relatively low power. 


86 







































































































Received QAM symbols prior to RS decoding 
16 QAM, BER2.5E-4 
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Figure 4.12. 16 QAM Signal Realization,= 6 kHz, R = 5 kbps, SNR = 7.4 dB 

BER = 2.5 E-4. 
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V. CONCLUSIONS 


The primary thrust of this thesis was to develop a computer simulation of an 
OFDM system for underwater acoustic communication using a PE based model of the 
acoustic channel. The simulation was validated by comparison of the results of the 
simulations with the known theoretical results established in this thesis and those from 
ocean experiments. The comparison of the results with what is expected from the theory 
was made in conjunction with presenting the results. Coatelan and Glavieux [Ref. 5] 
experimented with a OFDM communication system in shallow water acoustic channels 
with results as shown in Table 5.1. The QAM constellation size used by Coatelan and 
Glavieux [Ref. 1] is not known. Note that the water depth of the experiments is much 
shallower than for the acoustic channels used in the simulations of this thesis. The EEC 
code employed in the experiments, which uses rate I /2 convolutional coding with a 
constraint length of seven as well as a soft Viterbi decoder working on a 64 state trellis 
[Ref 1], is more advanced than that of the simulation, however the coding gain is likely 
within a couple of dB of the coding gain for the simulation as more than a few dB of 
coding gain improvement beyond 5 dB is difficult to obtain. Comparing the composite 
results presented earlier for the 8 and 16 QAM constellations and applying a margin of 6 
to 12 dB one can see that simulation results are comparable to the experimental results of 
Table 5.1. 


Trajismission 

distance 

Water 

depth 

SeaflooT 

type 

Emitted 

number 
of bit 

BER after 

demodulation 

BER after 
decoding 

estimated 

average 
level E^/Nq 

650 m 

25 m 

sludge 

18500 

3,7.10"^ 

< I.IO"'’ 

27,5 dB 

1920 m 

27 m 

sand 

18500 

2,4.10"2 

< 1.10*^ 

16 dB 

Table 1 . Multicarrier irammLision performances during two stations 


Table 5.1. Experimental Multicarrier Results [from Ref. 1 ]. 
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The benefits of the work are the verifieation of the feasibility of using OFDM as a 
method of underwater aeoustie communication and to provide a test-bed for future work 
to develop improved methods for the blocks that make up the OFDM system. 

There is significant opportunity for future work in this area. Many of the 
algorithms applied in the OFDM system are simple and could be improved to enhance the 
performance of the OFDM system. 

The modulation is DSBSC which uses twice the bandwidth of an Upper Side 
Band (USB) system. There are no obvious reasons why USB modulation could not be 
used to improve the spectral efficiency. This idea also has benefit in working on a multi¬ 
user application of OFDM. 

Communication system protocols must be established to make the system useful 
for real world application and appears to have been neglected so far in the area of 
underwater acoustic communication. The IEEE standard 802.11 and the European digital 
terrestrial TV broadcast uses OEDM and has established protocols that may be applicable 
to underwater communication systems. 

Synchronization algorithms have significant merit for future study and 
improvement. This study does not involve Doppler effects since it assumes that the 
transmitter and receiver are stationary. 

The application of the channel transfer function must be modified in modeling 
source or receiver motion. Erequency offset due to wave action and phase noise are not 
simulated either. Both of these are real concerns in development of an actual system and 
are worthy of simulation and analysis. The synchronization method must be modified in 
the simulation to account for phase issues. In addition the synchronization method is very 
computationally intensive and simplification of the system would benefit the simulation 
and any resulting real system. 

Channel estimation and equalization assumes a linear time invariant acoustic 
channel. The channel is known to be time varying with the rate of variance usually low 
but still of concern. Therefore a method of updating the channel estimation is necessary. 
Also while the estimation method used in the simulation is very robust there may be 
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better methods that use less bandwidth. FEC coding applied is elementary and can be 
improved. 

Thus far work to improve the system and the degree of reality of the simulation 
have been discussed. Of equal importance is the study of the key parameters on the 
system performance. The peak to average power ratio of the transmit signal is of concern 
in OFDM and has not been addressed in this thesis. 
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APPENDIX A. SIMULATION RESULTS 


The figures in this appendix represent the results of the simulation for all of the 
acoustic channels evaluated. The parameters of the system are provided in each figure as 
well as the main channel parameters. Each figure includes a zoomed inset which 
illustrates the BER performance near the target level of 10'"^. 



Figure A. 1. BER vs SNR for 8 QAM, fc = 8 kHz, R = 5 kbps. Range = 2 km. Rough 

= Om, Water Depth = 100m. 
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Figure A.2. BER vs SNR for 8 QAM, fc = 8 kHz, R = 5 kbps. Range = 4 km. Rough 

= Om, Water Depth = 100m. 



SNR (dB) 


Figure A.3. BER vs SNR for 8 QAM, fe = 10 kHz, R = 6 kbps. Range = 4 km. Rough 

= Om, Water Depth = 100m. 
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Figure A.4. BER vs SNR for 8 QAM, fc = 10 kHz, R = 6 kbps. Range = 2 km. Rough 

= Om, Water Depth = 100m. 



Figure A.5. BER vs SNR for 8 QAM, fe = 8 kHz, R = 5 kbps. Range = 2 km. Rough 

= Om, Water Depth = 100m. 
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Figure A.6. BER vs SNR for 8 QAM, fc = 12 kHz, R = 7.5 kbps. Range = 4 km. 

Rough = Om, Water Depth = 100m. 



Figure A.7. BER vs SNR for 16 QAM, fc = 6 kHz, R = 5 kbps. Range = 2 km. Rough 

= Om, Water Depth = 100m. 
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Figure A.8. BER vs SNR for 16 QAM, fc = 6 kHz, R = 5 kbps. Range = 2 km. Rough 

= 2m, Water Depth = 100m. 



Figure A..9 BER vs SNR for 16 QAM, fc = 6 kHz, R = 5 kbps. Range = 2 km. Rough 

= 4m, Water Depth = 100m. 
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Figure A. 10. BER vs SNR for 16 QAM, fc = 6 kHz, R = 5 kbps. Range = 4 km. Rough 

= Om, Water Depth = 100m. 


SNR = 2IdB, = 6kHz, R = 5kbps, Range = 4 km. Rough = 4, Water Depth = 100 m 



SNR (dB) 


Figure A. 11. BER vs SNR for 16 QAM, fe = 6 kHz, R = 5 kbps. Range = 4 km. Rough 

= 4m, Water Depth = 100m. 
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Figure A. 12. BER vs SNR for 16 QAM, fc = 6 kHz, R = 5 kbps. Range = 2 km. Rough 

= Om, Water Depth = 340m. 



Figure A. 13. BER vs SNR for 16 QAM, fc = 6 kHz, R = 5 kbps. Range = 4 km. Rough 

= 2m, Water Depth = 100m. 
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Figure A. 14. BER vs SNR for 16 QAM, fc = 6 kHz, R = 5 kbps. Range = 2 km. Rough 

= 2m, Water Depth = 340m. 


SNR = 21dB, f = 6kHz, R = 5kbps, Range = 4 km. Rough = 0, Water Depth = 340 m 



SNR (dB) 


Figure A. 15. BER vs SNR for 16 QAM, fc = 6 kHz, R = 5 kbps. Range = 4 km. Rough 

= Om, Water Depth = 340m. 
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Figure A. 16. BER vs SNR for 16 QAM, fc = 6 kHz, R = 5 kbps. Range = 4 km. Rough 

= 2m, Water Depth = 340m. 


SNR = 21dB, = 12kHz, R = 10kbps, Range = 2 km. Rough — 0, Water Depth = 100 m 



SNR(dB) 


Figure A. 17. BER vs SNR for 16 QAM, fc = 12 kHz, R = 10 kbps. Range = 2 km. 

Rough = Om, Water Depth = 100m. 
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SNR = 21 dB, f = 12kHz, R = 10kbps, Range = 4 km. Rough = 0, Water Depth = 100 m 



SNR (dB) 


Figure A.18. BER vs SNR for 16 QAM, fc = 12 kFIz, R = 10 kbps, Range = 4 km. 

Rough = Om, Water Depth = 100m. 



Figure A. 19. BER vs SNR for 16 QAM, fc = 8 kHz, R = 6.67 kbps. Range = 2 km. 

Rough = Om, Water Depth = 100m. 
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SNR = 21dB, = 8kHz, R = 6.6667kbps, Range = 4 km. Rough = 0, Water Deptli = 100 m 
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Figure A.20. BER vs SNR for 16 QAM, fc = 8 kHz, R = 6.67 kbps. Range = 4 km. 

Rough = Om, Water Depth = 100m. 



Figure A.21. BER vs SNR for 16 QAM, fe = 10 kHz, R = 8 kbps. Range = 2 km. 

Rough = Om, Water Depth = 100m. 
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Figure A.22. BER vs SNR for 16 QAM, fc = 10 kHz, R = 8 kbps. Range = 4 km. 

Rough = Om, Water Depth = 100m. 



Figure A.23. BER vs SNR for 32 QAM, fc = 10 kHz, R = 10 kbps. Range = 2 km. 

Rough = Om, Water Depth = 100m. 
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Figure A.24. BER vs SNR for 32 QAM, fc = 10 kHz, R = 10 kbps. Range = 4 km. 

Rough = Om, Water Depth = 100m. 



Figure A.25. BER vs SNR for 32 QAM, fe = 8 kHz, R = 8.33 kbps. Range = 2 km. 

Rough = Om, Water Depth = 100m. 
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Figure A.26. BER vs SNR for 32 QAM, fc = 8 kHz, R = 8.33 kbps. Range = 4 km. 

Rough = Om, Water Depth = 100m. 


SNR = 25dB, = 6kHz, R = 5kbps. Range = 2 km. Rough = 0, Water Depth = 100 m 



SNR (dB) 


Figure A.27. BER vs SNR for 32 QAM, fc = 6 kHz, R = 5 kbps. Range = 2 km. Rough 

= Om, Water Depth = 100m. 
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Figure A.28. BER vs SNR for 32 QAM, fc = 6 kHz, R = 5 kbps. Range = 2 km. Rough 

= 2m, Water Depth = 100m. 



Figure A.29. BER vs SNR for 32 QAM, fe = 6 kHz, R = 5 kbps. Range = 2 km. Rough 

= 4m, Water Depth = 100m. 
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Figure A.30. BER vs SNR for 32 QAM, fc = 6 kHz, R = 5 kbps. Range = 4 km. Rough 

= Om, Water Depth = 100m. 



Figure A.31. BER vs SNR for 32 QAM, fc = 6 kHz, R = 5 kbps. Range = 4 km. Rough 

= 2m, Water Depth = 100m. 
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Figure A.32. BER vs SNR for 32 QAM, fc = 6 kHz, R = 5 kbps. Range = 4 km. Rough 

= 4m, Water Depth = 100m. 



Figure A.33. BER vs SNR for 32 QAM, fe = 6 kHz, R = 5 kbps. Range = 2 km. Rough 

= Om, Water Depth = 340m. 
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SNR = 25dB, f = 6kHz, R = 5kbps, Range = 2 km. Rough = 2, Water Deptli = 340 m 



SNR (dB) 


Figure A.34. BER vs SNR for 32 QAM, fc = 6 kHz, R = 5 kbps. Range = 2 km. Rough 

= 2m, Water Depth = 340m. 



Figure A.35. BER vs SNR for 32 QAM, fc = 6 kHz, R = 5 kbps. Range = 4 km. Rough 

= Om, Water Depth = 340m. 
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Figure A.36. BER vs SNR for 32 QAM, fc = 6 kHz, R = 5 kbps. Range = 4 km. Rough 

= 2m, Water Depth = 340m. 



Figure A.37. BER vs SNR for 64 QAM, fc = 8 kHz, R = 10 kbps. Range = 2 km. 

Rough = Om, Water Depth = 100m. 
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Figure A.38. BER vs SNR for 64 QAM, fc = 8 kHz, R = 10 kbps. Range = 4 km. 

Rough = Om, Water Depth = 100m. 



Figure A.39. BER vs SNR for 64 QAM, fc = 6 kHz, R = 5 kbps. Range = 2 km. Rough 

= Om, Water Depth = 100m. 
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Figure A.40. BER vs SNR for 64 QAM, fc = 6 kHz, R = 5 kbps. Range = 2 km. Rough 

= Om, Water Depth = 100m. 



Figure A.41. BER vs SNR for 64 QAM, fc = 6 kHz, R = 5 kbps. Range = 2 km. Rough 

= 4m, Water Depth = 100m. 
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Figure A.42. BER vs SNR for 64 QAM, fc = 6 kHz, R = 5 kbps. Range = 4 km. Rough 

= Om, Water Depth = 100m. 



Figure A.43. BER vs SNR for 64 QAM, fc = 6 kHz, R = 5 kbps. Range = 4 km. Rough 

= 2m, Water Depth = 100m. 
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Figure A.44. BER vs SNR for 64 QAM, fc = 6 kHz, R = 5 kbps. Range = 4 km. Rough 

= 4m, Water Depth = 100m. 



Figure A.45. BER vs SNR for 64 QAM, fe = 6 kHz, R = 5 kbps. Range = 2 km. Rough 

= Om, Water Depth = 100m. 
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Figure A.46. BER vs SNR for 64 QAM, fc = 6 kHz, R = 5 kbps. Range = 2 km. Rough 

= 2m, Water Depth = 100m. 



Figure A.47. BER vs SNR for 64 QAM, fc = 6 kHz, R = 5 kbps. Range = 4 km. Rough 

= Om, Water Depth = 100m. 
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Figure A.48. BER vs SNR for 64 QAM, fc = 6 kHz, R = 5 kbps. Range = 4 km. Rough 

= 2m, Water Depth = 100m. 
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APPENDIX B. MMPE RESULTS 
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Figure B. 1 Transmission Loss, Roughness 0 m, Range 2 km 
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Figure B.2 Acoustic Paths, Roughness 0 m. Range 2 km,^^ 7 kHz 
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Figure B.3 Transmission Loss, Roughness 2 m, Range 2 km 
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Figure B.4 Aeoustic Paths, Roughness 2 m. Range 2 km,^^ 7 kHz 
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Figure B.5 Transmission Loss, Roughness 4 m, Range 2 km 



Figure B.6 Aeoustic Paths, Roughness 4 m. Range 2 km,^^ 7 kHz 
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Figure B.7 Transmission Loss, Roughness 0 m, Range 2 km 
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Figure B.8 Aeoustic Paths, Roughness 0 m. Range 2 km,/; 7 kHz 
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Figure B.9 Transmission Loss, Roughness 2 m, Range 2 km 
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Figure B. 10 Acoustic Paths, Roughness 0 m. Range 2 km,^^ 7 kFIz 
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Figure B. 11 Transmission Loss, Roughness 0 m, Range 2 km 



Figure B.12 Acoustic Paths, Roughness 0 m. Range 2 km,^^ 10 kHz 
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I'ransmission Loss (dB re Im) 



Figure B. 13 Transmission Loss, Roughness 0 m, Range 2 km 



Figure B.14 Acoustic Paths, Roughness 0 m. Range 2 km,^^ 12 kHz 
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Figure B. 15 Transmission Loss, Roughness 0 m, Range 2 km 
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Figure B.16 Acoustic Paths, Roughness 0 m. Range 2 km,^^ 14 kFlz 
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APPENDIX C. MATLAB CODE 


A, TRANSMITTER CODE 


1. ofdm sim xmitter.m 


%LT Tiger Pittman 
%UWA OFDM Thesis m-file 
%Last modified 4/16/01 
elear all, elose all 

fn = 0; %running figure number index 

synehdee = 0; 

%logie variable for synehronization, 0 == perfeet synehronization, 
%1 == synehronization using eyelie prefix eorrelator 



zero fraetion = .0250; %fraetion of subearriers that 

%oarry zeroes, near Fs/2 to avoid pseudo-aliasing effeets 

[N_bar, N_subo_aot, N_zero, guard_time, T, T_fft, Fs_fft, zero_fraotion, W, fe, 
m, SNR, delay spread] = parameters(zero_fraotion); 

%takes input from user to determine key parameters in the OFDM system 

%ineluding the number of aetive subearriers, guard time(eorrelates to eyelie 
prefix length) 

%OFDM symbol period(T, ineludes the FFT/IFFT interval and the guard time). 



noise = 1; 

%logie variable for noise addition 
noisefaetor = 0.35; 

%if noise == 1 then AWGN is added with following SNR to reeeived signal 
%otherwise no noise is added 


SNR_awgn = SNR*noise_faetor; 
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modulate = 2; 

%if modulate = 2 then double side band modulation is performed, otherwise 
%no modulation is performed 



decision = 2; 

%determines the information signal to be used 
%1 is for sinusoidal signals 
%2 is for a converted voice wave file 
%3 is for a ramp function 



%Signal generation 

pts = 2^13; %number of sample points (same as number of quant 

values) 

Fs=250; %samplmg frequency 

fl = 50; %tonal frequency 

f2 = 95; 

al = 1; %amplitude 

a2 = 0; 
factor = 20; 

[info sig, fin] = signal_generator(pts,Fs,fl,f2,al,a2,decision, factor, fin); 

%calls to get discrete time analog signal 



%Quantization 

bits = 2^3; %number of bits used in converting symbol(mteger) to 

binary 

disp([num2str(bits),' bit uniform quantization is performed resulting in ', 
num2str(2^bits),' quantization levels.']) 

[quant sig, codebook] = quantization(info_sig, bits); 

%performs UNIFORM quantization of info_sig to create a 

%discrete time, discrete valued signal, quant_sig(i) 
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%Conversion (quantized to binary) 

bit matrix = binary_conversion(quant_sig, bits); 

%pts by bits matrix of binary data 



%reshape (binary matrix to vector) 

bit_stream = reshape(bit_matrix.', 1, size(bit_matrix,l)*size(bit_matrix,2)); 
%converts bit matrix to bit stream for DMT use 

%note; bit matrix is transposed since the reshape command operates columnwise 
%but the data in bit matrix is oriented row wise 



%QAM and Reed Solomon parameters 

N = 2^m - 1; %Code word length for Reed Solomon symbol encoding 

q = m; %number of bits per QAM symbol, must be <= m 

Q = 2^q; %QAM constellation size, i.e. number of different 

QAM symbols 

%Note: as long as q is even, the QAM constellation is Gray coded (Desired) 
if(Q == 2^6); 

k = 47; %RS message length of Qary - QAM symbols 

elseif (Q == 2^5) 
k = 23; 

elseif (Q == 2M) 
k = 7; 

elseif (Q == 2^3) 
k = 3; 

elseif (Q == 2^2) 
k=l; 
else 

errorCQAM constellation size improperly determined') 
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%QAM Mapping (binary word to QAM index) 
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qam coded dec = qam_buffer(bit_stream, q); 

%buffers bit stream to q column matrix and converts 
%to decimal format vector with values 0 to Q-1 of length 
%ceil(length(bit_stream)/q) 



%Reed Solomon encoding 
if(Q >= 2^3) 

[rs code, zero symbols] = reedsolomon(qam_coded_dec, N, k); 
%encodes QAM symbols using Reed Solomon coding 
%rs_code is a matrix with each row of length N 
else 

rscode = qamcodeddec; 
zerosymbols = 0; 
end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%Data reshape (matrix to vector) 

code stream = reshape(rs_code, size(rs_code,l)*size(rs_code,2), 1); 
%reshapes rs code into column vector for QAM modulation 
%NOTE: the reshape command takes data columnwise. However rs_code 
%is not transposed in the reshape command. This effectively 
%interleaves the symbols since the subsymbols are not taken 
%in sequence 



%QAM Modulation 

[X, data carriers, zero fdl] = qam_modulator(code_stream, N bar, N_subc_act, 
N zerc, Q); 

%X is the N(=2*N_bar) row by data_blocks column matrix representing 
%the spectrum of the transmitted signal 
Fs_prime = 2*Fs_fft; 

%sampling rate doubles due to the doubling of the size of the data which 
%takes place when X is created conjugate symmetric 
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%IFFT (converts complex QAM symbols into real time domain sequence) 
X = inverse_fft(X); 

%N by data blocks matrix of real time domain samples 
%parallel output of the IFFT 



%Cyelie extension application 
L = round(guard_time*Fs_prime); 

%L = 2^ceil(log2(guard_time*Fs_prime)); %expected length of the 

impulse response of the chaimel 

%eomputed using the guard time and the sampling rate at the IFFT 

%this is the number of smaples in the guard time, or cyelie extension period, 
which must also be an integer 

x_cyclic = cye_extension(x, N_bar, L); 

%inserts a cyclic extension on each of the blocks(symbols) 



%Data reshape (2*N_bar by data_blocks matrix to vector) 
xmitsig = parlel_2_serial(x_cyclic); 

%serial time domain samples, length N*data blocks 



%Modulation (baseband signal modulation to carrier frequency) 
if (modulate == 1) 

[mod_xmit_sig, Fs_xmit, f_m, fn] = usb_modulator(xmit_sig, fc, Fs_prime, W, 
zero fraction, fn); 

%modulates baseband signal xmit sig to modulation frequency, fc 

%xmit_sig is interpolated to inerease Fs to Fs_xmitin order to modulate w/out 
aliasing 

elseif (modulate == 2) 

[mod_xmit_sig, Fs xmit, f_m, n_lpf2, len_xint, fn] = dsb_modulator(xmit_sig, fc, 
Fs_prime, W, zero fraction, fn); 

%modulates the baseband signal using double side band modulation 
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end 


2, parameters.m 

function [N bar, N subc act, N zerc, guard_time, T, T_fft, Fs_fft, zero_fraction, 
W, fc, m, SNR, delay spread] = parameters(zero_fraction) 

%This function uses input data to determine key parameters in the OFDM signal 
generation. 

%Last modified 7/3/01. 

R = 1000*(5.0+0/3);%input('What is the desired data rate (kbps)? '); 

Wmax = 1000*7; 
fc = 1000*6; 
delay_spread = le-3*25; 

Pe = .0001; 

%one dimensional symbol error probability, 2 D symbol error probability = 2*Pe 
SNR_gap3 = 20*logl0(sqrt(2)*erfmv(l-Pe)) 

SNR_gap = SNR_gap3 - 10*logl0(3) 

SNR = 21; 

b_Pe = log2(l + 10^(SNR/10)/10^(SNR_gap/10)); 

%bits per subsymbol that can be carried at given Pe, SNR gap, coding gain = 
margin dB and SNR 

guardtime = 4*delay_spread; 

%guard time required to prevent ISI (sec) 

T = 6*guard_time; 

%OFDM symbol period including the cyclic prefix 
Rs = 1/T; 

%OFDM symbol rate (symbols/sec) 
bofdm = R*T; 

%bits per OFDM symbol required to get desired rate 
Nsubcact = ceil(round(b_ofdm)/floor(b_Pe)); 

%number of subcarriers required for desired data rate at achievable b_Pe 
%these are data carrying subcarriers 
power2 = 2^ceil(log2(N_subc_act)); 

%determines the next power of 2 above N_subc_act 
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if (((power2-N_subc_act)/power2) >= zero_fraction) 

N_bar = power2; 
else 

N_bar = 2^ceil(log2(N_sube_act) + 1); 
end 

N_zerc = N_bar - N_sube_act; 

%number of zero subearriers 
zero_fraotion = N_zero/N_bar; 

%fraotion of subearriers that earry zeroes 
T fft = T - guard time; 

%FFT/IFFT period 
Fs_fft = N_bar/(T_fft); 

%tirst guess at FFT sampling rate in FFT interval 

%based on guard time = 4*expeeted delay time and an OFDM symbol period of 
%6* guard time 

%this is half the sampling rate required at the DAC and LPF sinee the 
%IFFT doubles the number of samples thereby effeetively interpolating by 2 
%and doubling the sampling rate 
n_samp = Fs_fft*T; 

%number of samples in the OFDM symbol period based on Fs fft 
%must also be an integer to ensure that the subcarriers are orthogonal 
if (isposint(n_samp)~=l) 
while (isposint(n_samp)~=l) 

Fs_fft = Fs_fft + 2; 
nsamp = Fs_fft*T; 
end 
end 

%revises the sampling rate slightly in the FFT interval such that there are an 
integer 

%number of samples in the OFDM symbol period 
T_fft = N_bar/Fs_fft; 

%revises the FFT/IFFT interval length based on the sampling rate 
%such that there will be an integer number of samples in the interval 
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guard time = T - T fft; 

%revises the guard time based on the new T fft and the known T 
format long e 

SC_del_f = (T-guardtime)^-1 
%subcarrier spaeing (Hz) 
format short e 

W_mmpe = SC_del_f''N_bar 
W = SC_del_PN_subo_aot; 

%bandwidth in Hz 
if (W > W max) 

disp('Unable to meet requirements with given bandwidth limitation') 
end 

%system parameters output 
m = floor(b_Pe); 
if (m==5|m==3) 

dispCQAM eonstellation is not Gray eoded') 
end 

disp('!!!!!!!!!!! RESULTS !!!!!!!!!!!!!!!') 

disp(['System bandwidth for DSBSC transmission is ',num2str(2*W/1000),' 
(kHz)']) 

fimax = fc + W; 

fimin = fe-W; 

if (fmax>7e3|fmin<200) 

disp('Frequency spectrum of system is outside the bandwidth of the known 
channel transfer function') 

end 

disp(['System Transmission Band for USBSC is ', num2str(fmin/1000),' (kHz) to ', 
num2str(fmax/1000),' (kHz).']) 

disp( ['Maximum Delay Spread is ', num2str(delay_spread*1000),' (msec)']) 

disp(['System uses ', num2str(N_subc_act),' active subcarriers with 
',num2str(N_bar),' total subcarriers.']) 

disp(['The subcarrier bandwidth is ', num2str(SC_del_f), ' (Hz).']) 
disp(['OFDM symbol duration is ', num2str(T*1000),' (msec).']) 
disp(['OFDM Symbol rate is ', num2str(Rs),' (symbols/sec).']) 
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disp( ['Number of bits per subcarrier is ',num2str(m),'.']) 

disp(['System bit rate is num2str(R/1000), ' (kbps).']) 

disp(['System Bandwidth Efficiency is ',nuni2str(R/(W*2)),'.']) 

disp(['System carries ', nuni2str(floor(b_Pe)*N_subc_act),' bits per OFDM 
symbol.']) 


3, signal generator.m 

function [x, fin] = signal_generator(pts,Fs,fl,f2,al,a2,decision, factor, fin) 
%This function generates a discrete time analog signal. 

%The signal is a vector of length pts, from the superposition of 
%2 unit amplitude sinusoidal signals of frequency fl and f2, sampled 
%at Fs Hz. 

%Fast modified 5/15/01 
%global fin 
if decision == 1 
n=0:pts-l; 

X = al *sin(2*pi*fl *n/Fs)+ a2*cos(2*pi*f2*n/Fs); %information signal 

spec_info_sig = fft(x,2048); 

fn = fin+l; figure(fn), plot(abs(spec_info_sig)) 

title('FFT of sampled information signal (info sig)') 

elseif decision ==2 

s = wavread('numnuts'); 

1 = round(length(s)/factor); 

x = s(l;l); 

elseif decision == 3 

X = [zeros(l,100), 0;.2;1, ones(l,100), l:-.2;0, zeros(l,100)]; 
spec_info_sig = fft(x,2048); 

fn = fin+l; figure(fn), plot(abs(spec_info_sig)) 
title('FFT of sampled information signal (info sig)') 
end 
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4, quantization.m 

function [symbol, codebook] = quantization(info_sig, bits) 

%Last modified 3/19/01 

M = 2^bits; %number of uniform quantization levels 
%must determine the method for partitioning to use in general 
maxsample = max(abs(info_sig)); 

%absolute maximum value of input signal 

part=-(M-l)/M*max_sample;2/M*max_sample:(M-l)/M*.9*max_sample; 

%partitions 

eodebook = -((M-l)*2 + l)/(2*M)*max_sample:2/M*max_sample:((M- 

l)*2+l)/(2*M)*max_sample; 

symbol = quantiz(info_sig, part); 

5, dec_2_bin,m 

funetion bit stream = dee_2_bin(qam_sym, bits, q) 

%Last modified 3/19/01 
bit_matrix = de2bi(qam_sym',q); 

%binary matrix length(qam_sym) by q 

bit_stream_long = reshape(bit_matrix, length(qam_sym)*q, 1); 

num zero bits = length(qam_sym)*q - floor((length(qam_sym)*q)/bits)*bits; 

%number of zero bits inserted in xmitter 

bit_stream = bit_stream_long(num_zero_bits+l ;length(bit_stream_long)); 

6, qam buffer.m 

function qam_dec = qam_buffer(bit_stream, q) 

%This function buffers the bit stream into a matrix form for QAM eneoding. 
%Last modified 3/19/01 

addlength = q*(ceil(length(bit_stream)/q))-length(bit_stream); 

%determines the number of zeroes needed to be added to bit stream 
%sueh that buffered matrix will be reetangular and inelude all of the data. 
bit_stream_long = [zeros(l,add_length),bit_stream]; 

%bit_stream with zeroes added to left end i.e. at 1,2,3... 
buffermat = reshape(bit_stream_long,length(bit_stream_long)/q,q); 
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%buffers bit stream into q column matrix for QAM modulation 
%data is reshaped columnwise 

qam_dec = bi2de(buffer_mat); %converts buffered binary data to 
corresponding 

%decimal integer where the integers are from 0 to Q-l(2^q - 1) 

%this reduces the sampling frequency since there are now less samples 
%least sig bit assumed in column 1 


7, reedsolomon.m 

function [code, zero symbols] = reedsolomon(qam_coded_dec, N, k) 

%This function performs Reed Solomon encoding of the decimal format 
%Gray coded Qary QAM symbols in qam coded dec 
%k is the message length in symbols, N is the code word length in symbols 
%Last modified 3/19/01 

full_blocks = floor(length(qam_coded_dec)/k); 

%number k length blocks filled with data 

zero symbols = (full_blocks + l)*k - length(qam_coded_dec); 

%number of zero blocks needed to fill final block 
qamcodeddeclong = [zeros(zero_symbols,l);qam_coded_dec]; 

%QAM code with zeroes at top of column vector 
qam_k_mat = reshape(qam_coded_dec_long,full_blocks+l,k); 

%full_blocks + 1 row by k column QAM symbol matrix for block coding 
code = encode(qam_k_mat, N, k, 'rs/decimal'); 

%full_blocks + 1 row by N column matrix of Qary-QAM encoded symbols 
%RS encoded 

8, qam modulator.m 

function [X, data carriers, zero_fill] = qam_modulator(code_stream, N bar, 
N_subc_act, N_zerc, Q) 

%Last modified 5/15/01 

[inphase, quad] = qaskenco(code_stream,Q); 

%encodes integer decimals into quadriture and inphase components of Q-ary 
QAM 
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[phase qam, mag qam] = cart2pol(inphase,quad); 

%converts QAM subsymbols to magnitude and phase(radians) for IFFT 
zero_carriers = N_zerc; %ceil(zero_fraction*N_bar); 

%number of zero subcarriers near Fs/2 
data_carriers = N_bar - zero carriers; 

%number of subcarriers with encoded data 
datablocks = ceil(length(phase_qam)/data_carriers); 

%number of blocks( i.e. symbols) neeeded for given N bar and zero fraction 
zero fill = data_blocks*data_carriers - length(phase_qam); 
colz = floor(zero_fill/data_blocks); 

%number of zeroes to be added into QAM symbol matrix for columns 2 ; 
datablocks 

col lz = zero fill - (data_blocks-l) * col z; 

%number of zeroes to be added to the first column in the QAM symbol matrix 

zero fill = [col lz, ones(l, data_blocks+2-l)*col_z]; 

phase cl = phase_qam(l:(data_carriers - col lz)); 

phase col = phase_qam(l+(data_carriers - col_lz):length(phase_qam)); 

col_phase_mat = reshape(phase_col, data_carriers - col_z, data_blocks-l); 

mag cl = mag_qam(l;(data_carriers - col lz)); 

mag col = mag_qam(l+(data_carriers - col_lz);length(mag_qam)); 

col_mag_mat = reshape(mag_col, data carriers - col z, data_blocks-l); 

tpcl = [phase_cl(l:col_lz).';zeros(l,col_lz)]; 

tpc2 = reshape(tpcl, 2*col_lz, 1); 

phaseclO = [tpc2; phase_cl(col_lz+l;length(phase_cl))]; 

tmcl = [mag_cl(l;col_lz).';zeros(l,col_lz)]; 

tmc2 = reshape(tmcl, 2*col_lz, 1); 

magclO = [tmc2; mag_cl(col_lz+l;length(mag_cl))]; 

for col = 1 : data blocks - 1 

tpcola = [col_phase_mat(l:col_z, col).'; zeros(l, col z)]; 
tpcolb = reshape(tpcola, 2*col_z, 1); 

phasecolO(:,col) = [tpcolb;col_phase_mat(col_z+l;size(col_phase_mat,l), col)]; 
tmcola= [col_mag_mat(l;col_z, col).'; zeros(l, col z)]; 
tmcolb = reshape(tmcola, 2*col_z, 1); 
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magcolO(:,col) = [tmcolb;col_mag_mat(col_z+l;size(col_mag_mat,l), col)]; 
end 

phase qam mat = [phaseclO, phasecolO]; 
mag qam mat = [mage 10, mageolO]; 
for k = 1 ;data_blocks 

X_half(: ,k) = zeros(N_bar+1,1); 
for 1=1 :data_carriers; 

X_half(i+l,k) = mag_qam_mat(i,k).*exp(j.*phase_qam_mat(i,k)); 
end 

X_half(l,k) = real(X_half(data_carriers+l,k)); 
X_half(data_carriers+l,k) = imag(X_half(data_carriers+l,k)); 

X(:,k) = [X_half(:,k); conj(flipud(X_half(2;N_bar,k)))]; 
end 

count = 0; 

for r = data_carriers+l :-l; 1 
for c = size(X_half,2):-l;l 
count = count + 1; 

X_pilot_half(count,l) = X_half(r,c); 
end 
end 

X_half(: ,data_bloeks+1) = zeros(N_bar+1,1); 

X_half(:,data_blocks+2) = zeros(N_bar+l,l); 

X_half(2;data_carriers+l,data_blocks+l) = X_pilot_half(l ; data carriers); 

X_half( 1 ,data_blocks+1) = real(X_half(data_carriers+1 ,data_bloeks+1)); 

X_half(data_carriers+1 ,data_bloeks+1) 
imag(X_half(data_carriers+1 ,data_blocks+1)); 

X_half(l ;data_carriers+l ,data_blocks+2) 
flipud(X_half( 1: data_carriers+1 ,data_bloeks+1)); 

for c = 1:2 

for r = 1 :data_carriers+l 
if X_half(r,2) == 0 
X_half(r,data_blocks+o) = 0; 
end 
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end 

end 

for k = data_blocks+l:data_blocks+2 
X(:,k) = [X_half(:,k); eonj‘(flipud(X_half(2;N_bar,k)))]; 
end 

X_pilots = X(:, data_blocks+l;size(X,2)); 

X = [X_pilots, X(:,l:data_blocks)]; 
for k = 1 :data_blooks+2 
for i=l;data_carriers; 
if(k <= datablocks) 

X_half(i+l,k) = mag_qam_mat(i,k).*exp(j.*phase_qam_mat(i,k)); 

%creates first half of complex symbol, other half is the 

conjugate 

%symmetric counterpart 
else 

X_half(i+l,k) = 1; 
end 
end 

X_half(l,k) = real(X_half(data_carriers+l,k)); 

%first channel is the real part of the last symbol 
if k <= datablocks 

X_half(data_carriers+l,k) = imag(X_half(data_carriers+l,k)); 
end 

X(:,k) = [X_half(:,k); conj(fiipud(X_half(2;N_bar,k)))]; 
end 


9, inverse fft.m 

function x_n = inverse_fft(X) 

%Last modified 4/13/01 
x_n = ifft(X); 

%N time domain samples, real valued from symmetric fireq spectrum, X 
%x_n is a matrix the same size as X, representing the time domain samples 
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%to be transmitted for each block 
max_imag = max(max(imag(x_n))); 
if (max_imag>= 1 e-14) 

error('Time domain output from IFFT not real to machine accuracy') 
end 


10, cyc_extension,m 

function x_cycex = cyc_extension(x, N_bar,L) 
%Last modified 3/19/01 
for i = 1 ;size(x,2) 

x_cycex(:,i) = [x(2*N_bar-L+l;2*N_bar,i);x(:,i)]; 
end 


11, parlel_2_serial,m 

function x_ser = parlel_2_serial(x_cyclic) 

%This function converts the parallel output of the IFFT to serial form. 
%x_ser(l) represents the last sample to be transmitted 
%Last modified 3/19/01 

X = reshape(x_cyclic, 1, size(x_cyclic,l)*size(x_cyclic,2)); 
x_ser = real(x); 

%real command necessary to correct for small imaginary 
%components in x from machine error 

%converts symbol from parallel to serial w/ x(N,data_blocks)leading the signal 
%and x(l,l) trailing 


12, dsb_modulator,m 

function [mod_sig, Fs mod, f_m, n_lpf2, len_xint, fn] = dsb_modulator(xmit_sig, 
fc, Fs_prime, W, zero fraction, fn); 

%This function modulates the baseband signal, xmit sig to the carrier frequency, 
fc, using 

%double sideband suppressed carrier (DSBSC) modulation. 
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%The baseband signal is interpolated as necessary to meet Nyquist frequency 
requirements 

%Last modified 4/5/01 


%plots 

N=8192; %fft length 

f_p=(0:N-l)*Fs_prime/N; %index for plotting fffs 



%ERROR CHECK TO ENSURE NO AEIASING 


if (fc < W) 

error(Aliasing will occur due to improper fc and W combination, fc must be 
greater than W.') 

end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%%%%%%%%%%%%% 


% UPSAMPEING TO INCREASE THE SAMPEING RATE 

K= 1.0; 

% margin above the Nyquist rate 
I = ceil(2*K*(fc+W)/Es_prime); 

%interpolation constant 

[x int, fin] = interp_mod(xmit_sig, I, fin); 

%interpolated xmit signal, baseband, DSB 
lenxint = length(x_int); 

Esmod = I*Es_prime; 

% sampling frequency of xmitted and modulated signal 
f_m = (0;N-l)*Es_mod/N; 

%frequency index for interpolated signal spectrum 



%POST INTERPOEATION EPE 
fracpass2 = 1.2; 

% fraction of W for pass band 
fracstop2 = 1.375; 

% fraction of W for stop band 
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digf_pass2 = 2*pi*fracpass2*W/(Fs_mod); % relative pass frequency 
digf_stop2 = 2*pi*fracstop2*W/(Fs_mod); % relative stop frequency 
xtion2 = digf_stop2-digf_pass2; %transition region width 

fpass2 = fracpass2*W/(.5*Fs_mod); 

%fraction of Fs/2 for cutoff frequency in FIR filter from firl.m 

n lpfZ = round( 8*pi/xtion2 ); %hamming filter order 

if (isposint(n_lpf2/2) ~= 1) %ensures fdter has odd number of 

coefficients 

n_lpf2 = n_lpf2+l; 

end 

[Nh2,Dh2] = firl( n_lpf2 , fpass2 ); %Hamming LP filter numerator and 

denominator coefficients 

x_int_lpf = conv(Nh2, x_int); 



mod = cos(2*pi*fc/Fs_mod*[l:length(x_int_lpf)]); 
modsig = mod.*x_int_lpf; 

mod_sig_spec = fft(mod_sig, N);%, 2^ceil(log2(length(mod_sig)))); 
xmit_sig_spec = fft(xmit_sig, N);%, 2^ceil(log2(length(xmit_sig)))); 


B, RECEIVER CODE 


1, ofdm sim receiver.m 

%LT Tiger Pittman 
%OFDM Simulation; receiver 
%Last Modified 5/16/01 
depth = 0; 

%depth = 0 corresponds to shallow water model, 100m depth 
%depth = 1 corresponds to deep water model, 340m depth 
range = 1; 

%range 0 corresponds to 2 km transmission range, 1 to 4 km range 
rough = 0; 
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%can be 0,2 or 4 

%0 for all allowed fc, 2 for deep water and fc = 7 only, 
%2 and 4 for shallow water at fc = 7 only 
recvdsig = mod_xmit_sig; 



%ADDITIVE WHITE GAUSSIAN NOISE 
if noise == I 

%Ps=(sum(recvd_sig( 1; 1000) A2)/1000); 

Ps=(sum(recvd_sig( 1 ;length(recvd_sig)) A2)/length(recvd_sig)); 

Pn = Ps/((I0^(SNR_awgn/I0))); 

awgn = randn(I,length(recvd_sig))*sqrt(Pn); 

recvd_sig = recvd_sig + awgn; 

% verification of SNR in time domain 
Pnoise=(sum(awgn( 1:1000) A2)/1000); 
SNR_est=10*logl0(Ps/Pnoise); % checks SNR 
display(SNR_est) 
end 



%DEMODULATION, (includes low pass filtering and decimation) 
if (modulate == 1) 

[recvd sig, fin] = rcvr_usb_demod(recvd_sig, Fs xmit, fc, mod_xmit_sig, 
W, Fs_prime, fin); 

elseif (modulate == 2) 

[recvd sig, fin] = rcvr_dsb_demod(recvd_sig, Fs xmit, fc, mod_xmit_sig, W, 
Fs_prime, n lpfZ, len xint, synchdec, fin); 

end 



%performs time and frequency synchronization 
blocklength = size(x_cyclic,l); 
num_blocks = size(X,2); 
if synchdec == 1 
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[sync_sig, synch, fn] = synchronizer(recvd_sig, block length, num_blocks, 
N_bar, L, SNR_awgn, fn); 

else 

sync_sig = recvd_sig; 
end 



%Serial to parallel conversion 

par_sig_r = ser_2_parlel(sync_sig, L, N_bar); 

%parallel blocks(symbols) of received data with cyclic extension 
%still attached at top of each symbol(column) 



%Cyclic extension removal 

par = par_sig_r(L + l:2*N_bar + L, :); 



%FFT of received signal 
X_r_eb = fft(par); 

%FFT of received symbols, includes channel equalization blocks 



%Channel application 

[X_r_ebch, rcvr_depth, fn] = mmpe_channel(X_r_eb, N_bar, N_subc_act, fc, 
depth, range, rough, fn); 



%performs channel equalization 

X_r = equalization(X_r_ebch, X, N_subc_act, N_bar, fn); 



%QAM Demodulator 

[code, inphase r, quad_r, fn] = qam_demodulator(X_r, N bar, zero_fraction, Q, 
data_carriers, N, zero_fdl, fn); 

%demodulates X_r to give Reed Solomon encoded integers, 0;Q-1 
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codedeint = deinterleaver(code, N); 
%deinterleaves the received symbols 



if(Q >= 2^3) 

qam_sym = rsdecoder(code_deint, N, k, zero_symbols); 
%decodes RS encoded symbols 

else 

qamsym = codedeint; 
end 



bit_stream_r = dec_2_bin(qam_sym, bits, q); 

%reshapes qam sym to q column matrix for conversion to binary 
%removes zero bits inserted in xmitter 



buffer_r = reshape(bit_stream_r, bits, length(bit_stream_r)/bits)'; 
%reshapes bit_stream_r to bits column matrix for conversion to 
%decimal integers for signal reconstruction 



quantint = bi2de(buffer_r); 

%converts binary words of length bits to decimal integer 



infosigr = codebook(quant_int + 1); 



signals2(info_sig, info_sig_r, code, code stream, Q, X_r, fin) 
BER = error_check(bit_stream, bit_stream_r); 


146 


2, rcvr dsb demodulator.m 


function [dec_f_rsig, fn] = rcvr_dsb_demod(recvd_sig, Fs_xmit, fc, 
mod xmit sig, W, Fs_prime, n_lpf2, len_xint, synchdec, fn); 

%This function demodulates the DSBSC signal from the ehannel transfer funetion 

%Last modified 5/9/01 



%DEMODULATION 

demodsig = eos(2*pi*fe*(l;length(reevd_sig))/Fs_xmit); 
%demodulating signal to eombine upper side band speetrums to reereate 
%baseband signal 

demodreesig = reovd_sig.*demodsig; 

%demodulated time domain signal 



N = 2*8192; %length(demod_reo_sig); 
f = (0:N-l)*Fs_xmit/N; %index for plotting fft's 

%frequenoy index of reeeived signal speetrum 



%LOW PASS FILTERING 


fraepass = 1.075; 

% fraetion of W for pass band 
fraestop = 1.2; 

% fraetion of W for stop band 

digf_pass=2*pi*fraopass*W/(Es_xmit); % relative pass frequeney 

digf_stop=2*pi*fraestop*W/(Es_xmit); % relative stop frequeney 

if (digf_stop>=.99*pi) 

error('Stop band exeeeds .99 of the Sampling frequeney') 
end 

xtion=digf_stop-digf_pass; %transition region width 

fpass = fraepass* W/(.5 *Es_xmit); 

%fraetion of Es/2 for eutoff frequeney in EIR filter from firl.m 
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n lpf = round( 8*pi/xtion ); %hamming filter order 

if (isposint(n_lpf/2) ~= 1) %ensures filter has odd number of 

coefficients 

nlpf = n_lpf+l; 

end 

[Nh,Dh] = firl( n lpf, fpass); %Hamming LP filter numerator and 

denominator coefficients 

fdt_rsig = conv(Nh, demod_rec_sig); 

if synchdec == 0 

timedelaymod = nlpfZ; 

timedelaydemod = nlpf; 

time_delay_tot = fioor(time_delay_mod/2) + fioor(time_delay_demod/2); 

filt_rsig = filt_rsig(l, time_delay_tot+l;len_xint+time_delay_tot); 

%TH1S COMMAND PERFORMS PERFECT SYNCHRONIZATION 

end 



%DEC1MAT10N 


D = round(Fs_xmit/Fs_prime); 
if (isposint(D)~=l) 

error('Decimation factor not an integer') 
end 

if (isposint(length(filt_rsig)/D)~= 1) 
ifD==2 

filt_rsig = [filt_rsig, 0]; %appends zero to signal so that it will decimate 
evenly 

elseif D >= 3 

lenl = length(filt_rsig) + 1; 
len2 = length(filt_rsig) + 2; 
len3 = length(filt_rsig) + 3; 
if isposint(lenl/D) 
filt rsig = [filt rsig, 0]; 
elseif isposint(len2/D) 
filt rsig = [fdt_rsig, 0 0]; 
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else 

filt rsig = [filt rsig, 0 0 0]; 
end 
end 
end 

ifD>=5 

error('Interpolation/Decimation factors greater than 4 requires code adjustment 
in rcvr dsb demod.m') 

end 

dec_mat = reshape(fdt_rsig, D, length(fdt_rsig)/D); 
dec_f_rsig = dec_mat(l,:); 

%decimated received signal, returns the sampling frequency to Fs_prime = 
Fs_xmit/D 

f_p = (0;N-l)*Fs_prime/N; %mdex for plotting fffs 

%frequency index of received signal spectrum 


3, synchronizer,!!! 


function [syncsig, synch, fin] = synchronizer(recvd_sig, blocklength, 
num_blocks, N_bar, L, SNR_awgn, fn) 

%Correlator synchronizer 

%last modified 6/21/2001 



%ADDITIVE WHITE GAUSSIAN NOISE 


%Ps=(sum(recvd_sig( 1: 1000) A2)/1000); 

Ps=(sum(recvd_sig( I ;Iength(recvd_sig)) A2)/length(recvd_sig)); 
%n=0:I;N-I; 

%S=sin(2*pi*16.6*n/PS)+sin(2V*I7.9''n/PS); 

Pn = Ps/((I0^(SNR_awgn/I0))); 

awgn = randn(I,length(recvd_sig))*sqrt(Pn); 

% verification of SNR in time domain 
Pnoise=(sum(awgn( I; 1000) A2)/1000); 
SNR_est_synch=I0*IogI0(Ps/Pnoise); % checks SNR 
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display(SNR_est_synch) 
lengnoisel = 1000; 
leng_noise2 = 1500; 

r = [sqrt(Pn)*randn(l, leng noisel), reevd_sig, sqrt(Pn)*randn(l, leng_noise2) ]; 



for k = 1 ;length(r)-2*N_bar-2*L+l %loops over entire length of data 

r_sig = r(l, k:k+L-l); 

r_del_oonj_sig = conj(r(l, k+2*N_bar:k+2*N_bar+L-l)); 
corr_sig(k) = max(xoorr(r_sig, r_del_conj_sig)); 
end 

[synch, fin] = peakfmder(corr_sig, fin); 

%last block identifies the ofdm block length 

%previous blocks identify the discrete time prior to the beginning of the data 
blocks 

%i.e. synch(l) is the beginning of the first data block and synch(end) is the ofdm 
block length 

%performance check 

block_length_error = abs(block_length - synch(length(synch))); 
peak count error = num blocks - (length(synch) - 1); 
block start error = abs(leng_noisel + 1 - synch(l)); 

disp(['Block length error = ', num2str(block_length_error),' discrete time units.']) 

disp(['Peaks missed = ', num2str(peak_count_error),'.']) 

sync_sig = r(synch(l)+l : synch(l) + 1 + num_blocks*block_length - 1); 


a. peakfinder.m 

function [synch, fin] = peakfmder(corr_sig, fin) 

%last modified 5/16/01 
peakindex = 0; 

%peaks = find(corr_sig>mn + 2.5 *sd); 
fracnum = 0; 

tl = .4; t2 = .6; dt = .05; %threshholds and differential 

for fraction = tl ;dt:t2 
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fracnum = fracnum + 1; 

%fraction = 0.7; 

peaks = find(corr_sig > fraction*max(corr_sig)); 
count = 0; 
peaknum = 0; 

while count < length(peaks) 
same_peak = 1; 

peaknum = peaknum + 1; %index for number of peaks found 
lookstart = count + 1; 
while same_peak == 1; 
count = count + 1; 
if (count >= length(peaks)) 
break 
end 

if (peaks(l+count) - peaks(count) == 1) 
same_peak = 1; 
else 

same_peak = 0; 
end 
end 

%while and if loop determine the indices of corr sig 
corresponding to the 1st peak which is where 

%the first OFDM symbol begins 

%count-l after while loop is the last index of the 1st peak, 
nextpeakstart = count + 1; 

[mag_peak(peaknum), peak(peaknum)] = 

max(corr_sig(peaks(lookstart);peaks(count))); 

peak_index(fracnum, peaknum) = peak(peaknum) + 

peaks(lookstart) - 1; 

end 

end 

%disp(fraction) 

%disp(peak_index) 
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numtruepeaks = 0; 
for 1 = 1 :size(peak_index,2) 
if peak_index(l,l) ~= 0 

matches = find(peak_index(l,l) == 

peak_index(2:size(peak_index,l),:)); 

if length(matches) == size(peak_index,l) - 1 

numtruepeaks = numtruepeaks + 1; 

truepeakindex(numtruepeaks) = peak_index(l,l); 

end 

end 

end 

%break 

block length = truepeakindex(2) - truepeakindex(l) - 1; 
synch = [truepeakindex, block length]; 

fn=fn+l; figure(fn), plot(corr_sig), title('Correlator Synchronizer Output'), 
hold on 

plot(truepeakindex, corr_sig(truepeakindex), 'r*') 
for k = tl;dt:t2 

line([l, length(corr_sig)], [k*max(corr_sig), k*max(corr_sig)]) 
end 

hold off 

4, mmpe channeLm 

function [X_ch, rcvr depth, fin] = mmpe_channel(X_r_eb, N_bar, N_subc_act, fc, 
depth, range, rough, fn) 

%This function applies the transfer function from the MMPE model to the 
received 

%QAM symbols. X_r_ebch(:,i) = X_r_eb(:,i).*H(z), where i is the OFDM block 
index 

% 

%Last modified 6/4/01 
ifN_bar== 1024 
if fc ~= 6e3 

error('Carrier frequency does not correspond to model carrier frequency') 
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end 

if (depth == 0 & range == 0 & rough == 0) 
load f7r02k.mat; 

elseif (depth == 0 & range == 0 & rough == 2) 
load f7r22k.mat; 

elseif (depth == 0 & range == 0 & rough == 4) 
load f7r42k.mat; 

elseif (depth == 0 & range == 1 & rough == 0) 
load f7r04k.mat; 

elseif (depth == 0 & range == 1 & rough == 2) 
load f7r24k.mat; 

elseif (depth == 0 & range == 1 & rough == 4) 
load f7r44k.mat; 

elseif (depth == 1 & range == 0 & rough == 0) 
load f7r02kd.mat; 

elseif (depth == 1 & range == 0 & rough == 2) 
load f7r22kd.mat; 

elseif (depth == 1 & range == 1 & rough == 0) 
load f7r04kd.mat; 

elseif (depth == 1 & range == 1 & rough == 2) 
load f7r24kd.mat; 

end 

elseif N bar == 2048 


if (fc == 8e3 & range == 0) 
load fl0r02k.mat; 
elseif (fc == 8e3 & range == 1) 
load fl0r04k.mat; 
elseif (fc == 10e3 & range == 0) 
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load fl2r02k.mat; 
elseif (fc == 10e3 & range == 1) 
load fl2r04k.mat; 
elseif (fc == 12e3 & range == 0) 
load fl4r02k.mat; 
elseif (fc == 12e3 & range == 1) 
load fl4r04k.mat; 
end 
end 

r = size(press,l); %number of depths in model, first half are in 

acoustic channel, i.e. water 

[v dep_index] = max(sum(abs(press(l:r/2,:)),2)); %dep_index indetifies 

the depth with the greatest 

%return over all frequencies 

H = press(dep_index,:); %places the receiver at the optimal depth 

%H = press(10,:); 

H_half = zeros(length(H)+l,l); 
for k = 1 :length(H) 

H_half(k+1) = H(k); 
end 

H_half(l) = real(H_half(length(H_half))); 

H_half(length(H_half)) = imag(H_half(length(H_half))); 

H full = [H half; conj(fhpud(H_half(2;N_bar)))]; 

%transfer function to be applied to received QAM symbols 
h = flipud(ifft(H_full)); 

rcvr_depth = (l-depth)*(dep_index/(r/2))*100 + depth*(dep_index/(r/2))*340; 
%receiver depth in meters 
for k = 1 :size(X_r_eb,2) 

X_ch(:,k) = X_r_eb(:,k).*H_full; 

End 

5. equalization.m 

function X_r = equalization(X_r_eb, X, N subc act, N_bar, fin) 
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%performs channel equalization using pilot blocks 
%last modified 5/15/01 

H_est_mat_real = X_r_eb(l:N_subc_act+l, 1:2)./ ... 

X(l:N_subc_act+l, 1:2); 

H_est_mat_imag = X_r_eb(size(X_r_eb,l) - N_subc_act + 1: size(X_r_eb,l), 

1 : 2 )./... 

X(size(X,l) - N_subc_act + 1: size(X,l), 1:2); 

%estimate of channel freqency response from 2 unit magnitude OFDM blocks 

%estimate is the average of the 2 columns and is used for channel equalization 

%real is estimate for first N_subc_act subcarriers, and imag is for last N_subc_act 
subcarriers 

%equalization is only performed over the active subcarriers 
H_est_real = mean(H_est_mat_real.').'; 

H_est_imag = mean(H_est_mat_imag.').'; 

%2*N_subc_act column vector representing the channel estimate over the actice 
subcarriers 

X_r = X_r_eb(:, 3:size(X_r_eb,2)); 

%removes channel estimation blocks assuming the first 2 blocks are the pilot 
blocks and no others 

forkl = l:size(X_r,2); 

X_r_real(:,kl) = X_r(l:N_subc_act+l,kl)./H_est_real; 

X_r_imag(:,kl) = X_r(size(X_r,l) - N_subc_act + 1: 

size(X_r, l),kl)./H_est_imag; 

X_r(:,kl) = [X_r_real(:,kl); zeros(2*N_bar - 2*N_subc_act - 1,1); 

X_r_imag(:,kl)]; 

end 

%equalizes the data for the channel 


6, qam demodulator.m 

function [code, inphase_r, quad_r, fn] = qam_demodulator(X_r, N bar, 
zero_fraction, Q, data_carriers, N, zero fill, fn) 

%Last modified 5/11/01 

X_half_r = X_r(l:data_carriers+1,:); 

%strips conjugate symmetric and zero subcarriers from data 
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for k = 1; size(X_half_r,2) 
for i = l;data_carriers 
mag_r(i,k) = abs (X_half_r(i+l,k)); 

%determines the magnitude of the QAM symbol 
phase_r(i,k) = angle(X_half_r(i+l,k)); 

%determmes the phase of the QAM symbol 
end 

mag_r(data_earriers,k) = sqrt(X_half_r(l,k)^2 + 

real(X_half_r(data_earriers+1 ,k))^2); 

%eombmes the Re(N_bar+l) and the Imag(N_bar+l) that were separated 

%in the qam modulator to get the magnitude if the final symbol 

phase_r(data_oarriers,k) = atan2(real(X_half_r(data_oarriers+1 ,k)), 
X_half_r(l,k)); 

%oombines the Re and Im at N bar + 1 to get the phase for the final 

symbol 

%real eommand neeessary to aeeount for small imaginary eomponent 
resulting from 

%maohine error 

end 

%result of these 2 loops is a data earrier row matrix representing the QAM 
modulated 

%symbols where eaeh eolumn is the kth symbol and eaeh row is the ith 
subsymbol 



tmagl(:,l) = mag_r(l;2:zero_fill(l)*2,l); 
tphasel(:,l) = phase_r(l;2:zero_fdl(l)*2,l); 
for eol = 2;size(X_r,2) 

tmagool(:,ool-l) = mag_r(l;2;zero_fill(ool)*2,ool); 
tphaseool(:,ool-l) = phase_r(l:2:zero_fdl(ool)*2,ool); 
end 

mag_r_veel(:,l) = [tmagl(;,l); mag_r(zero_fdl(l)*2+l;length(mag_r),l)]; 
phase_r_veel(:,l) = [tphasel(:,l); phase_r(zero_fill(l)*2+l;length(phase_r),l)]; 
for eol = 2:size(X_r,2) 
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[tmagcol(:,col-l); 


mag_matcol(:,col-1) = 

mag_r(zero_lill(col)*2+l;length(mag_r),col)]; 

phase_matcol(:,col-1) = [tphasecol(:,col-l); 

phase_r(zero_fill(col)*2+1 ;length(phase_r),col)]; 

end 

mag_r_veccol = reshape(mag_matcol, size(mag_matcol,l)*size(mag_matcol,2), 

1 ); 

phaserveccol = reshape(phase_matcol, 

size(phase_matcol,l)*size(phase_matcol,2), 1); 

mag_r_vec_short = [mag_r_vec 1; mag r veccol]; 

phase r vec short = [phase_r_vecl; phase_r_veccol]; 



[inphase_r, quad_r] = pol2cart(phase_r_vec_short, mag_r_vec_short); 
%converts the mag and phase values to in phase and quadriture components 
%for QAM demodulation 

minmax = [min(inphase_r), max(inphase_r); min(quad_r), max(quad_r)]; 
code = qaskdecomod(inphase_r, quad r, Q, minmax); 

%decodes the QAM symbols into integers 0:Q-1 


a. qaskdecomod.m 

function [msg]=qaskdecomod(x,y,m,minmax) 

%QASKDECO Decodes QASK mapped signal back to integer symbols. 

% MSG = QASKDECO(X, Y, M) decodes the message signal MSG 
from the 

% in-phase component X and quadrature component Y with the M-ary 

% number M. M must equal 2^K with K being an positive integer. The 

% minimum and the maximum numbers are assumed to be the same as 
in the 

% output from QASKENCO. 

% 

% MSG = QASKDECO(X, Y, M, MINMAX) decodes the information 
where 

% the maximum and minimum values of X and Y are given in 

MINMAX 
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% 

by the form: 


% 

X min 

X max 

% 

MINMAX = 1 


% 

Y min 

Y max 

% 



% 

See also QASKENCO. 

% 

East modified 5/11/01 


scale = 0; %logic statement to seale reeeived symbols, 1 = yes, o/w = 
no 

%if seale = 0 

% dispCNO SCALING OF RECEIVED QAM CONSTEEEATION 
PEREORMED IN QAM DECODER') 

%end 

error(nargehk( 1,4,nargin)); 
if (nargin <= 2) | (nargin > 4) 

error('Wrong number of input arguments.') 
end; 

M = m; 
if isstr(m) 

error(['QASK eonstellation has ehanged to use Gray eode. ', ... 
sprintf('\n'),... 

'Use COMMUPDTC'EIEENAME") to update your SIMUEINK 
model'... 

sprintf('\n'),... 

'before the simulation.']); 

end; 

K = log(m) ./ log(2); 
ifm~=2^K 

error('M must equal to 2^K where K is an integer.') 
end; 

XX = constlay(K, 1); 

[leny, lenx] = size(xx); 
if (nargin == 4) 
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[xmat,ymat] = qaskenco(M); 

%xmat and ymat are the inphase and quad eomponents, respeetively, of 
the Mary QAM 

%eonstellation 

MAT = [xmat.';ymat.']; 

%2 by M matrix of veetors in 2D representing the M Qam symbols in 
the eonstelation 

for k=l;length(x) 

E = [x(k); y(k)]*ones(l,length(xmat)) - MAT; 

[v, index] = min([l,l]*abs(E)); 
msg(k) = index-1; 
end 
end; 

[x_m, x_n] = size(x); 

[y_m, y_n] = size(y); 
if min(x_m, x_n) * min(y_m, y_n) 1 
error('Input X and Y must be vectors.'); 
end; 

if x_m > x_n 
msg = msg'; 
end; 

7, deinterleaver.m 

function deint_mat = deinterleaver(code, N) 

%East modified 3/19/01 
rows = length(code)/N; 

%number of rows in matrix for RS decoding 
deint_mat = reshape(code, rows, N); 

%matrix of symbols deinterleaved, ready for RS decoder 
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8, rsdecoder.m 


function decode = rsdecoder(code_deint, N, k, zero symbols) 

%Last modified 3/19/01 

deco = rsdeco(code_deint, N, k, 'rs/decimal'); 

%message length k columns matrix 

deco2 = reshape(deco, 1, size(deco,l)*size(deco,2)); 

%converts to row vector 

decode = deco2(zero_symbols+l;length(deco2)); 

%removes zeroes added in transmitter 

%zero symbols passed from xmitter function reedsolomon.m 


9, error_check.m 

function ber = error_check(bs, bsr) 

%This function determines the Bit Error Rate (BER) for the system. 
%last modified 7/2/01 
be = 0; 

for 1 = 1 :length(bs) 
be = be + (bs(l)~=bsr(l)); 
end 

ber = be/1; 


C. PERFORMANCE ANALYSIS CODE 

1, ofdm sim control func.m 

function [BER, SNR est] = ofdm_sim_control_func(SNR, fc, R, rough, range, 
depth) 

%Control m-file for ofdm simulation 
fn = 0; 

modulate = 2; %DSBSC modulation logic control value 

noise = 1; 
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for avn =1:10 
count = 0; 

for noisefactor = .15;.01;.65 
count = count + 1; 

synchdec = 0; %synchronization logic control value, 0 == NO, 1 

== YES 

[zero_fill, Q, N, data_carriers, N_bar, N_subc_act, zero_symbols, k, bits, 
info_sig, code_stream, bit_stream, X, mod_xmit_sig, ... 

Fs_prime, Fs xmit, fn, block_length, zero_fraction, SNR_awgn, chan_app, 
W, n lpfZ, len_xint, F, q, codebook] = ... 

ofdm_sim_xmitter_funo(synchdec, noise, noise factor, fn, R, fc, SNR); 

[SNR_est(count,avn), BER(count,avn)] = ofdm_sim_receiver_func(W, 
zero_fdl, Q, N, data_carriers, N_bar, N_subc_act, zero_symbols, k, bits, info_sig, 

codestream, bit_stream, X, mod_xmit_sig, Fs_prime, Fs_xmit, fn, 
block_length, zero fraction, SNR_awgn, ... 

depth, range, rough, chan_app, modulate, noise, fc, synchdec, n_lpf2, 
len xint, F, q, codebook); 

end 

end 


2, ofdm_performance,m 

%OFDM performance analysis m-file 

%runs ofdm sim control func m-file over all mmpe channels 

%last modified 6/28/01 


warning off 
clear all 



%%% 


%8-ary QAM PERFORMANCE 

SNR= 17 

fc = 8e3 

R=5e3 

count = 0; 
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for range = 0:1 %0 == 2 km, 1 == 4 km 

count = count + 1; 
rough = 0; 
depth = 0; 

[BERw, SNR estw] = ofdm_sim_control_func(SNR, fc, R, rough, range, 
depth); 

BER(:,:,count) = BERw; 

SNR_est(:,:,count) = SNR_estw; 
end 

clear BERw SNR_estw 
fc = 10e3 
R = 6e3 

for range = 0:1 %0 == 2 km, 1 == 4 km 

count = count + 1; 
rough = 0; 
depth = 0; 

[BERw, SNR_estw] = ofdm_sim_control_func(SNR, fc, R, rough, range, 
depth); 

BER(:,:,count) = BERw; 

SNR_est(:,:,count) = SNR_estw; 
end 

clear BERw SNR_estw 
fc = 12e3 
R=7.5e3 

for range = 0:1 %0 == 2 km, 1 == 4 km 

count = count + 1; 
rough = 0; 
depth = 0; 

[BERw, SNR_estw] = ofdm_sim_control_func(SNR, fc, R, rough, range, 
depth); 

BER(:,:,count) = BERw; 

SNR_est(:,:,count) = SNR_estw; 
end 
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clear BERw SNR estw 



%%% 

% 16-ary QAM PERFORMANCE 
SNR = 21 
fc = 6e3 
R = 5e3 
depth = 0; 

for range = 0:1 %0 == 2 km, 1 == 4 km 

for rough = 0:2:4 
count = count + 1; 

[BERw, SNR_estw] = ofdm_sim_control_func(SNR, fc, R, rough, range, 
depth); 

BER(:,:,count) = BERw; 

SNR_est(:,:,count) = SNR_estw; 
end %rough 
end %range 
clear BERw SNR_estw 
depth = 1; 

for range = 0:1 %0 == 2 km, 1 == 4 km 

for rough = 0:2:2 
count = count + 1; 

[BERw, SNR_estw] = ofdm_sim_control_func(SNR, fc, R, rough, range, 
depth); 

BER(:,:,count) = BERw; 

SNR_est(:,:,count) = SNR_estw; 
end %rough 
end %range 
clear BERw SNR_estw 
fc = 12e3 
R= 10e3 

for range = 0:1 %0 == 2 km, 1 == 4 km 

count = count + 1; 
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rough = 0; 
depth = 0; 

[BERw, SNR estw] = ofdm_sim_control_func(SNR, fc, R, rough, range, 
depth); 

BER(:,:,eount) = BERw; 

SNR_est(:,:,eount) = SNR_estw; 
end 

elear BERw SNR_estw 
fc = 8e3 

R= (6+2/3)* le3 

for range = 0:1 %0 == 2 km, 1 == 4 km 

count = count + 1; 
rough = 0; 
depth = 0; 

[BERw, SNR estw] = ofdm_sim_control_func(SNR, fc, R, rough, range, 
depth); 

BER(:,:,count) = BERw; 

SNR_est(;,:,count) = SNR_estw; 
end 

clear BERw SNR_estw 
fc = 10e3 
R = (8+0/3)* le3 

for range = 0:1 %0 == 2 km, 1 == 4 km 

count = count + 1; 
rough = 0; 
depth = 0; 

[BERw, SNR_estw] = ofdm_sim_control_func(SNR, fc, R, rough, range, 
depth); 

BER(:,:,count) = BERw; 

SNR_est(:,:,count) = SNR_estw; 
end 

clear BERw SNR estw 
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%%% 

%32-ary QAM PERFORMANCE 
SNR = 25 
fc = 10e3 
R=(10+0/3)Me3 

for range = 0:1 %0 == 2 km, 1 == 4 km 

count = count + 1; 
rough = 0; 
depth = 0; 

[BERw, SNR_estw] = ofdm_sim_control_func(SNR, fc, R, rough, range, 
depth); 

BER(:,:,count) = BERw; 

SNR_est(:,:,count) = SNR_estw; 
end 

clear BERw SNR_estw 
fc = 8e3 

R=(8+l/3)Me3 

for range = 0:1 %0 == 2 km, 1 == 4 km 

count = count + 1; 
rough = 0; 
depth = 0; 

[BERw, SNR_estw] = ofdm_sim_control_func(SNR, fc, R, rough, range, 
depth); 

BER(:,:,count) = BERw; 

SNR_est(:,;,count) = SNR_estw; 
end 

clear BERw SNR_estw 
fc = 6e3 

R= (5+0/3)* le3 
depth = 0; 

for range = 0:1 %0 == 2 km, 1 == 4 km 

for rough = 0:2:4 
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count = count + 1; 

[BERw, SNR_estw] = ofdm_sim_control_func(SNR, fc, R, rough, range, 
depth); 

BER(;,:,count) = BERw; 

SNR_est(:,;,count) = SNR_estw; 
end %rough 
end %range 
clear BERw SNR_estw 
fc = 6e3 

R= (5+0/3)* le3 
depth = 1; 

for range = 0:1 %0 == 2 km, 1 == 4 km 

for rough = 0:2:2 
count = count + 1; 

[BERw, SNR estw] = ofdm_sim_control_func(SNR, fc, R, rough, range, 
depth); 

BER(:,:,count) = BERw; 

SNR_est(:,:,count) = SNR_estw; 
end %rough 
end %range 
clear BERw SNR_estw 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 



%64ary QAM PERFORMANCE 
SNR = 27 
fc = 8e3 


R= (10+0/3)* le3 

for range = 0:1 %0 == 2 km, 1 == 4 km 

count = count + 1; 
rough = 0; 
depth = 0; 

[BERw, SNR_estw] = ofdm_sim_control_func(SNR, fc, R, rough, range, 
depth); 
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BER(:,:,count) = BERw; 

SNR_est(;,;,count) = SNR_estw; 
end 

clear BERw SNR_estw 
fe = 6e3 

R= (5+0/3)* le3 
depth = 0; 

for range = 0:1 %0 == 2 km, 1 == 4 km 

for rough = 0:2:4 
eount = eount + 1; 

[BERw, SNR estw] = ofdm_sim_oontrol_func(SNR, fc, R, rough, range, 
depth); 

BER(:,:,count) = BERw; 

SNR_est(:,:,count) = SNR_estw; 
end %rough 
end %range 
clear BERw SNR_estw 
fe = 6e3 

R= (5+0/3)* le3 
depth = 1; 

for range = 0:1 %0 == 2 km, 1 == 4 km 

for rough = 0:2:2 
count = count + 1; 

[BERw, SNR estw] = ofdm_sim_eontrol_fune(SNR, fc, R, rough, range, 
depth); 

BER(:,:,count) = BERw; 

SNR_est(:,:,count) = SNR_estw; 
end %rough 
end %range 
clear BERw SNR_estw 
warning on 

save perform BER SNR_est 
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3, perform plotter.m 


%m-file to plot BER vs SNR est for performance results 
clear all 
load perform 
fn = 0; 


close all 


%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%8 QAM 
for 1 = 1 ;6 

mn_ber8(:,l) = mean(BER(:,:,l),2); 

mn_snr8(:,l) = mean(SNR_est(:,:,l),2); 
end 
fnc = 0 
fnc = fnc+1 
figure(fnc) 

plot(mn_snr8,mn_ber8) 
title('8 QAM Composite Results') 
xlabel('SNR (dB)'), ylabel(’BER') 
legendCr, '2', '3', '4', '5', '6') 
%%%%%%%%%%%%%%%%%%%%%%%%% 


%16 QAM 


for 1 = 7:22 

mn_berl6(:,l) = mean(BER(:,:,l),2); 
mn_snrl6(:,l) = mean(SNR_est(:,:,l),2); 
end 

fnc = fnc+1 
figure(fnc) 

plot(mn_snr 16,mn_ber 16) 
title('16 QAM Composite Results') 
xlabel('SNR (dB)'), ylabel('BER') 


legendCr, '2', '3', '4', '5', '6', '7', '8', '9', '10', 'll', '12', '13', '14', '15', '16') 


%%%%%%%%%%%%%%%%%%%%%%%%%%% 
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%32 QAM 


for 1 = 23:36 

mn_ber32(:,l) = mean(BER(:,:,l),2); 
mn_snr32(:,l) = mean(SNR_est(:,:,l),2); 
end 

fnc = fnc+1 
figure(fnc) 

plot(mn_snr32,mn_ber32) 
title('32 QAM Composite Results') 
xlabel(’SNR (dB)'), ylabel(’BER') 


legendCr, '2', '3', '4', '5', '6', '7, '8', '9', '10', 'll', '12', '13', '14') 


%%%%%%%%%%%%%%%%%%%%%%%%% 


%64 QAM 


for 1 = 37:48 

mn_ber64(:,l) = mean(BER(:,:,l),2); 
mn_snr64(:,l) = mean(SNR_est(:,:,l),2); 
end 

fnc = fnc+1 
figure(fnc) 

plot(mn_snr64,mn_ber64) 
title('64 QAM Composite Results') 
xlabel('SNR (dB)'), ylabel('BER') 


legendCr, '2', '3', '4', '5', '6', '7, '8', '9', '10', 'll' 


' 12 ') 


%break 



%%% 


%8-ary QAM PERFORMANCE 
SNR= 17 
fc = 8e3 
R=5e3 


%%% 

fn = fin + 1; ligure(fn), orient tall 
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plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R 

',num2str(R/l 000), 'kbps, Range = 2 km. Rough = 0, Water Depth =100 m']) 

xlabel('SNR (dB)'), ylabel('BER') 

%%% 

fn = fn + 1; tigure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R 

',num2str(R/1000),'kbps, Range = 4 km. Rough = 0, Water Depth =100 m']) 

xlabel('SNR (dB)'), ylabel('BER') 

%%% 

fc = 10e3 

R = 6e3 

%%% 

fn = fn + 1; ligure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R 
',num2str(R/1000),'kbps, Range = 2 km. Rough = 0, Water Depth =100 m']) 

xlabel('SNR (dB)'), ylabel('BER') 

%%% 

fn = fn + 1; ligure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
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hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R = 
',num2str(R/l 000), 'kbps, Range = 4 km. Rough = 0, Water Depth =100 m']) 

xlabel(’SNR (dB)'), ylabel('BER') 

%%% 

fc = 12e3 

R = 7.5e3 

%%% 

fn = fn + 1; tigure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R = 

',num2str(R/1000),'kbps, Range = 2 km. Rough = 0, Water Depth =100 m']) 

xlabel('SNR (dB)'), ylabel('BER') 

%%% 

fn = fn + 1; ligure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R = 

',num2str(R/1000),'kbps, Range = 4 km. Rough = 0, Water Depth =100 m']) 

xlabel('SNR (dB)'), ylabel('BER') 

%%% 



%%% 


% 16-ary QAM PEREORMANCE 
SNR = 21 
fc = 6e3 
R=5e3 
%%% 
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fn = fn + 1; figure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R 

',num2str(R/l 000), 'kbps, Range = 2 km. Rough = 0, Water Depth =100 m']) 

xlabel(’SNR (dB)'), ylabel('BER') 

%%% 

fn = fn + 1; ligure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R 

',num2str(R/1000),'kbps, Range = 2 km. Rough = 2, Water Depth =100 m']) 

xlabel('SNR (dB)'), ylabel('BER') 

%%% 

fn = fn + 1; ligure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R 

',num2str(R/1000),'kbps, Range = 2 km. Rough = 4, Water Depth =100 m']) 

xlabel('SNR (dB)'), ylabel('BER') 

%%% 

fn = fn + 1; ligure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,:,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R 

',num2str(R/1000),'kbps, Range = 4 km. Rough = 0, Water Depth =100 m']) 
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xlabel('SNR (dB)'), ylabel(’BER') 

%%% 

fn = fn + 1; figure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R = 

',num2str(R/l 000), 'kbps, Range = 4 km. Rough = 2, Water Depth =100 m']) 

xlabel('SNR (dB)'), ylabel('BER') 

%%% 

fn = fn + 1; ligure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fo/1000),'kHz, R = 

',num2str(R/1000),'kbps, Range = 4 km. Rough = 4, Water Depth =100 m']) 

xlabel('SNR (dB)'), ylabel('BER') 

%%% 

%depth = 0; 

%for range = 0:1 %0 == 2 km, 1 == 4 km 

% for rough = 0:2:4 
% eount = eount + 1; 

% [BERw, SNR_estw] = ofdm_sim_oontrol_funo(SNR, fe, R, rough, range, 
depth); 

% BER(:,:,count) = BERw; 

% SNR_est(:,:,count) = SNR_estw; 

%end %rough 
%end %range 
%clear BERw SNR_estw 
%%% 

fn = fn + 1; ligure(fn), orient tall 
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plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R 

',num2str(R/l 000), 'kbps, Range = 2 km. Rough = 0, Water Depth = 340 m']) 

xlabel('SNR (dB)'), ylabel('BER') 

%%% 

fn = fn + 1; tigure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R 

',num2str(R/1000),'kbps, Range = 2 km. Rough = 2, Water Depth = 340 m']) 

xlabel('SNR (dB)'), ylabel('BER') 

%%% 

fn = fn + 1; tigure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R 

',num2str(R/1000),'kbps, Range = 4 km. Rough = 0, Water Depth = 340 m']) 

xlabel('SNR (dB)'), ylabel('BER') 

%%% 

fn = fn + 1; tigure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R 

',num2str(R/1000),'kbps, Range = 4 km. Rough = 2, Water Depth = 340 m']) 

xlabel('SNR (dB)'), ylabel('BER') 
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%%% 
fc = 12e3 
R= 10e3 
%%% 

fn = fn + 1; figure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), niean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R 

',nuni2str(R/l 000), 'kbps, Range = 2 km. Rough = 0, Water Depth =100 m']) 

xlabel(’SNR (dB)’), ylabel(’BER') 

%%% 

fn = fn + 1; ligure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R 

',num2str(R/1000),'kbps, Range = 4 km. Rough = 0, Water Depth =100 m']) 

xlabel('SNR (dB)'), ylabel(’BER') 

%%% 

fc = 8e3 

R = (6+2/3)* le3 
%%% 

fn = fn + 1; ligure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R 
',num2str(R/1000),'kbps, Range = 2 km. Rough = 0, Water Depth =100 m']) 

xlabel(’SNR (dB)'), ylabel(’BER') 
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%%% 

fn = fn + 1; figure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R = 
',num2str(R/l 000), 'kbps, Range = 4 km. Rough = 0, Water Depth =100 m']) 

xlabel(’SNR (dB)'), ylabel('BER') 

%%% 

fc = 10e3 

R= (8+0/3)* le3 

%%% 

fn = fn + 1; ligure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R = 

',num2str(R/1000),'kbps, Range = 2 km. Rough = 0, Water Depth =100 m']) 

xlabel('SNR (dB)'), ylabel('BER') 

%%% 

fn = fn + 1; ligure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R = 

',num2str(R/1000),'kbps, Range = 4 km. Rough = 0, Water Depth =100 m']) 

xlabel('SNR (dB)'), ylabel('BER') 

%%% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
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%32-ary QAM PERFORMANCE 
SNR = 25 
fc = 10e3 
R= (10+0/3)* le3 
%%% 

fn = fn + 1; figure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R 

',num2str(R/l 000), 'kbps, Range = 2 km. Rough = 0, Water Depth =100 m']) 

xlabel(’SNR (dB)’), ylabel(’BER') 

%%% 

fn = fn + 1; ligure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_o = ',num2str(fo/1000),'kHz, R 

',num2str(R/1000),'kbps, Range = 4 km. Rough = 0, Water Depth =100 m']) 

xlabel(’SNR (dB)'), ylabel(’BER') 

%%% 

fe = 8e3 

R=(8+l/3)*le3 

%%% 

fn = fn + 1; ligure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_o = ',num2str(fo/1000),'kHz, R 
',num2str(R/1000),'kbps, Range = 2 km. Rough = 0, Water Depth =100 m']) 
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xlabel('SNR (dB)'), ylabel(’BER') 

%%% 

fn = fn + 1; figure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R 
',num2str(R/l 000), 'kbps, Range = 4 km. Rough = 0, Water Depth =100 m']) 

xlabel('SNR (dB)'), ylabel('BER') 

%%% 

fc = 6e3 

R= (5+0/3)* le3 

%%% 

fn = fn + 1; ligure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_o = ',num2str(fc/1000),'kHz, R 

',num2str(R/1000),'kbps, Range = 2 km. Rough = 0, Water Depth =100 m']) 

xlabel('SNR (dB)'), ylabel('BER') 

%%% 

fn = fn + 1; ligure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_o = ',num2str(fo/1000),'kHz, R 

',num2str(R/1000),'kbps, Range = 2 km. Rough = 2, Water Depth =100 m']) 

xlabel('SNR (dB)'), ylabel('BER') 

%%% 

fn = fn + 1; ligure(fn), orient tall 
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plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R 

',num2str(R/l 000), 'kbps, Range = 2 km. Rough = 4, Water Depth =100 m']) 

xlabel('SNR (dB)'), ylabel('BER') 

%%% 

fn = fn + 1; tigure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R 

',num2str(R/1000),'kbps, Range = 4 km. Rough = 0, Water Depth =100 m']) 

xlabel('SNR (dB)'), ylabel('BER') 

%%% 

fn = fn + 1; ligure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R 

',num2str(R/1000),'kbps, Range = 4 km. Rough = 2, Water Depth =100 m']) 

xlabel('SNR (dB)'), ylabel('BER') 

%%% 

fn = fn + 1; ligure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R 

',num2str(R/1000),'kbps, Range = 4 km. Rough = 4, Water Depth =100 m']) 

xlabel('SNR (dB)'), ylabel('BER') 
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%%% 
fc = 6e3 

R= (5+0/3)* le3 
%%% 

fn = fn + 1; figure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), niean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R 
',nuni2str(R/l 000), 'kbps, Range = 2 km. Rough = 0, Water Depth = 340 m']) 

xlabel(’SNR (dB)’), ylabel(’BER') 

%%% 

fn = fn + 1; ligure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R 

',num2str(R/1000),'kbps, Range = 2 km. Rough = 2, Water Depth = 340 m']) 

xlabel('SNR (dB)'), ylabel(’BER') 

%%% 

fn = fn + 1; ligure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R 

',num2str(R/1000),'kbps, Range = 4 km. Rough = 0, Water Depth = 340 m']) 

xlabel('SNR (dB)'), ylabel(’BER') 

%%% 

fn = fn + 1; ligure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
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hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R = 
',num2str(R/l 000), 'kbps, Range = 4 km. Rough = 2, Water Depth = 340 m']) 

xlabel('SNR (dB)'), ylabel('BER') 

%%% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 



%64ary QAM PERFORMANCE 
SNR = 27 
fc = 8e3 

R = (10+0/3)Me3 
%%% 

fn = fin + 1; tigure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R = 

',num2str(R/1000),'kbps, Range = 2 km. Rough = 0, Water Depth =100 m']) 

xlabel('SNR (dB)'), ylabel('BER') 

%%% 

fn = fn + 1; tigure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R = 

',num2str(R/1000),'kbps, Range = 4 km. Rough = 0, Water Depth =100 m']) 

xlabel('SNR (dB)'), ylabel('BER') 

%%% 

fc = 6e3 
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R= (5+0/3)* le3 

%%% 

fn = fn + 1; figure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R 
',num2str(R/l 000), 'kbps, Range = 2 km. Rough = 0, Water Depth =100 m']) 

xlabel('SNR (dB)'), ylabel('BER') 

%%% 

fn = fn + 1; figure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fo/1000),'kHz, R 

',num2str(R/1000),'kbps, Range = 2 km. Rough = 2, Water Depth =100 m']) 

xlabel('SNR (dB)'), ylabel('BER') 

%%% 

fn = fn + 1; ligure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_o = ',num2str(fo/1000),'kHz, R 

',num2str(R/1000),'kbps, Range = 2 km. Rough = 4, Water Depth =100 m']) 

xlabel('SNR (dB)'), ylabel('BER') 

%%% 

fn = fn + 1; ligure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
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hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R 
',num2str(R/l 000), 'kbps, Range = 4 km. Rough = 0, Water Depth =100 m']) 

xlabel(’SNR (dB)'), ylabel('BER') 

%%% 

fn = fn + 1; tigure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R 

',num2str(R/1000),'kbps, Range = 4 km. Rough = 2, Water Depth =100 m']) 

xlabel('SNR (dB)'), ylabel('BER') 

%%% 

fn = fn + 1; ligure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R 

',num2str(R/1000),'kbps, Range = 4 km. Rough = 4, Water Depth =100 m']) 

xlabel('SNR (dB)'), ylabel('BER') 

%%% 

fc = 6e3 

R = (5+0/3)* le3 
%%% 

fn = fn + 1; ligure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R 
',num2str(R/1000),'kbps, Range = 2 km. Rough = 0, Water Depth = 340 m']) 

xlabel('SNR (dB)'), ylabel('BER') 
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%%% 

fn = fn + 1; figure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R 
',num2str(R/l 000), 'kbps, Range = 2 km. Rough = 2, Water Depth = 340 m']) 

xlabel(’SNR (dB)'), ylabel('BER') 

%%% 

fn = fn + 1; ligure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(:,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R 

',num2str(R/1000),'kbps, Range = 4 km. Rough = 0, Water Depth = 340 m']) 

xlabel('SNR (dB)'), ylabel('BER') 

%%% 

fn = fn + 1; ligure(fn), orient tall 
plot(SNR_est(:,:,fn), BER(;,:,fn), 'b.') 
hold on 

plot(mean(SNR_est(:,;,fn),2), mean(BER(:,:,fn),2), 'r-') 
hold off 

title(['SNR = ',num2str(SNR), 'dB, f_c = ',num2str(fc/1000),'kHz, R 

',num2str(R/1000),'kbps, Range = 4 km. Rough = 2, Water Depth = 340 m']) 

xlabel('SNR (dB)'), ylabel('BER') 
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